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1  Introduction 


Modern  space  surveillance  requires  fast  and  accurate  orbit  predictions  for  myriads  of  objects 
in  a  broad  range  of  Earth  orbits.  Conventional  Special  Perturbations  orbit  propagators, 
based  on  numerical  integration  of  the  osculating  equations  of  motion,  are  accurate  but  ex¬ 
tremely  slow  (typically  requiring  100  or  more  steps  per  satellite  revolution  to  give  good 
predictions).  Conventional  General  Perturbations  orbit  propagators,  based  on  fully  analyt¬ 
ical  orbit  theories  like  those  of  Brouwer,  are  faster  but  contain  large  errors  due  to  inherent 
approximations  in  the  theories.  New  orbit  generators  based  on  Semianalytic  Satellite  Theory 
(SST)  have  been  able  to  approach  the  accuracy  of  Special  Perturbations  propagators  and 
the  speed  of  General  Perturbations  propagators. 

SST  has  been  originated  by  P.  J.  Cefola  and  his  colleagues,  whose  names  are  in  the  refer¬ 
ences  at  the  end  of  this  document.  The  theory  is  scattered  throughout  the  listed  conference 
preprints,  published  papers,  technical  reports,  and  private  communications.  Our  purpose  in 
this  document  is  to  simplify,  assemble,  unify,  and  extend  the  theory.  This  document  includes 
truncation  algorithms  and  corrects  misprints  in  our  earlier  work  [Danielson,  Neta,  and  Early, 
1994], 

SST  represents  the  orbital  state  of  a  satellite  with  an  equinoctial  element  set  (oq, . . . ,  a^). 
The  first  five  elements  oq, . . .  ,05  are  slowly  varying  in  time.  The  sixth  element  06  is  the 
mean  longitude  A  and  so  is  rapidly  varying. 

SST  decomposes  the  osculating  elements  di  into  mean  elements  a*  plus  a  small  remainder 
which  is  27t  periodic  in  the  fast  variable: 

di  Oj  T  ,  •  •  •  ,  06,  t)  (1) 

(Here  we  use  hats  to  distinguish  the  elements  of  the  osculating  ellipse  from  the  elements  of 
the  averaging  procedure.  The  values  of  a  free  index  are  assumed  to  be  obvious  from  the 
context;  e.  g.  ,  here  i  can  have  the  values  1,  2,  3,  4,  5,  or  6,  so  (1)  represents  6  equations.) 
The  mean  elements  oq  are  governed  by  ordinary  differential  equations  of  the  form 

^  =  nSi 6  +  Ai(ai, . . . ,  a5,  t )  (2) 

Here  t  is  the  time,  n  is  the  (mean)  mean  motion,  and  8^  is  the  Kronecker  delta  (i.  e.  , 
<5 16  =  826  =  83Q  =  $46  =  $56  =  0,  $66  =  !)•  The  short-periodic  variations  rg  are  expressable  in 
Fourier  series  of  the  form 

OO 

Vi  =  (Q1>  •  •  •  >  «5,  t)  cos  jA  +  Si  (eq,  ...,a5,t)  sin  jX]  (3) 

j= 1 

Having  formulas  for  the  mean  element  rates  H,,  we  can  integrate  the  mean  equations  (2) 
numerically  using  large  step  sizes  (typically  1  day  in  length).  The  formulas  for  the  Fourier 
coefficients  Cj  and  Sj  in  (3)  also  only  need  to  be  evaluated  at  the  integrator  step  times. 
Values  of  the  osculating  elements  di  at  request  times  not  coinciding  with  the  integrator  step 
times  can  be  computed  from  (1)  using  interpolation  formulas. 

In  subsequent  chapters  we  will  outline  the  methods  of  derivation  and  give  explicit  formulas 
for  the  terms  Ai,  Cj ,  Sj  corresponding  to  various  perturbing  forces. 


3 


2  Mathematical  Preliminaries 


2.1  Equinoctial  Elements 

The  generalized  method  of  averaging  can  be  applied  to  a  wide  variety  of  orbit  element 
sets.  The  equinoctial  elements  were  chosen  for  SST  because  the  variational  equations  for 
the  equinoctial  elements  are  nonsingular  for  all  orbits  for  which  the  generalized  method  of 
averaging  is  appropriate-namely,  all  elliptical  orbits. 

In  this  chapter  we  give  an  overview  of  the  equinoctial  elements,  which  are  osculating 
(even  though  they  do  not  have  hats).  They  are  discussed  in  more  detail  in  [Broucke  and 
Cefola,  1972],  [Cefola,  Long,  and  Holloway,  1974],  [Long,  McClain,  and  Cefola,  1975],  [Cefola 
and  Broucke,  1975],  [McClain,  1977  and  1978],  and  [Shaver,  1980]. 


2.1.1  Definition  of  the  Equinoctial  Elements 

There  are  six  elements  in  the  equinoctial  element  set: 


CL\ 

&2 

«3 

04 

0,5 

06 


a  semimajor  axis 

components  of  the  eccentricity  vector 

components  of  the  ascending  node  vector 
A  mean  longitude 


h 

k 

P 

q 


The  semimajor  axis  a  is  the  same  as  the  Keplerian  semimajor  axis.  The  eccentricity  vector 
has  a  magnitude  equal  to  the  eccentricity  and  it  points  from  the  central  body  to  perigee. 
Elements  h  and  k  are  the  g  and  f  components,  respectively,  of  the  eccentricity  vector  in 
the  equinoctial  reference  frame  defined  below.  The  ascending  node  vector  has  a  magnitude 
which  depends  on  the  inclination  and  it  points  from  the  central  body  to  the  ascending  node. 
Elements  p  and  q  are  the  g  and  f  components,  respectively,  of  the  ascending  node  vector  in 
the  equinoctial  reference  frame. 

There  are  actually  two  equinoctial  element  sets:  the  direct  set  and  the  retrograde  set.  As 
the  names  imply,  the  direct  set  is  more  appropriate  for  direct  satellites  and  the  retrograde  set 
is  more  appropriate  for  retrograde  satellites,  ft  is  possible,  however,  to  use  direct  elements 
for  retrograde  satellites  and  vice  versa,  and  for  non-equatorial  satellites  this  presents  no 
problem.  For  equatorial  satellites  there  are  singularities  which  must  be  avoided  by  choosing 
the  appropriate  equinoctial  element  set.  For  direct  elements 


lim 


(1) 


while  for  retrograde  elements 


(2) 


For  each  equinoctial  element  set  there  are  three  associated  vectors  (f,  g,  w)  which  define 
the  equinoctial  reference  frame.  These  vectors  form  a  right-handed  orthonormal  triad  with 
the  following  properties: 
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Figure  1:  Direct  Equinoctial  Reference  Frame 


Figure  2:  Retrograde  Equinoctial  Reference  Frame 


1.  f  and  g  lie  in  the  satellite  orbit  plane. 

2.  w  is  parallel  to  the  angular  momentum  vector  of  the  satellite. 

3.  The  angle  between  f  and  the  ascending  node  is  equal  to  the  longitude  of  the  ascending 
node. 

This  leaves  two  choices  for  f  and  g,  one  associated  with  the  direct  element  set  and  one 
associated  with  the  retrograde  element  set.  The  two  sets  of  (f,  g,  w)  are  illustrated  in 
Figures  1  and  2.  In  the  Figures,  and  throughout  this  document,  ( x,y,z )  denote  a  set  of 
Cartesian  coordinates  whose  origin  moves  with  the  center  of  mass  of  the  central  body  and 
whose  axes  are  nonrotating  with  respect  to  inertial  space. 


2.1.2  Conversion  from  Keplerian  Elements  to  Equinoctial  Elements 

If  a,  e,  i,  hi,  u,  M  denote  the  conventional  Keplerian  element  set  then  the  equinoctial  ele¬ 
ments  are  given  by: 


a 

h 

k 

P 

q 

X 


e  sin  (a;  +  IQ) 
e  cos  (a;  +  /hi) 


tan 


tan 


(i) 

(i) 


sin  hi 
cos  hi 


M  +  ou  +  m 


(1) 


The  quantity  /  is  called  the  retrograde  factor  and  has  two  possible  values: 


I 


+1  for  the  direct  equinoctial  elements 
— 1  for  the  retrograde  equinoctial  elements 


(2) 


There  are  two  auxiliary  longitudes  associated  with  the  equinoctial  element  set:  the  ec¬ 
centric  longitude  F  and  the  true  longitude  L.  They  are  related  to  the  Keplerian  eccentric 
anomaly  E  and  true  anomaly  /  by  the  equations: 


F  =  E  +  u  +  m  (3) 

L  =  f  +  uj  +  m  (4) 


These  auxiliary  longitudes  are  used  in  converting  from  equinoctial  elements  to  position  and 
velocity.  In  addition,  certain  perturbations  are  modeled  with  Fourier  series  expansions  in  F 
or  L. 
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2.1.3  Conversion  from  Equinoctial  Elements  to  Keplerian  Elements 

In  order  to  convert  from  equinoctial  to  Keplerian  elements,  it  is  first  necessary  to  compute 
an  auxiliary  angle  d,  which  is  defined  by: 


sin  C  = 


'h 2  +  k2 


cos  d  = 

The  Keplerian  elements  are  then  given  by: 


' h 2  +  k2 


=  VWT¥ 
n-i\ 


^ ^  +  2 1  arctan  i Jp2  +  q2 


sin  12  = 


cos  12  = 


P 

/ p 2  +  q 2 

q 


LU  =  d  -  /ft 

M  =  A  -  d 

where  /  is  defined  by  (2. 1.2-2). 

The  Keplerian  eccentric  and  true  anomalies  are  given  by: 

E  =  F-C 

f  =  l- d 


2.1.4  Conversion  from  Equinoctial  Elements  to  Position  and  Velocity 

The  first  step  in  converting  from  equinoctial  elements  to  position  and  velocity  is  to  determine 
the  equinoctial  reference  frame  basis  vectors  (f,  g,  w).  Their  components  in  the  (x,  y,  z) 
system  are 


1  1  -  p2  +  q2 

1  +  P  + 1  [  -2  Ip 

T  2  Ipq 

T— T— -a  (1  +p2  -  q2)I 

!+?-  +  «-[  2  q 

1  '  2  p 

1 +p2  +  q2  _  (i_p2^g2)/ 


1  +p2  +  q2 


7 


The  second  step  is  to  find  the  eccentric  and  true  longitudes  F  and  L,  respectively.  To 
find  the  eccentric  longitude  F,  one  must  solve  the  equinoctial  form  of  Kepler’s  Equation  (see 
Section  7.1)): 

A  =  F  +  h  cosF  —  ksinF  (2) 

Then  define  two  auxiliary  quantities,  the  Kepler  mean  motion  n  and  a  quantity  called  b : 


n 

b 


1  +  VI  -  h2  -  k2 


(3) 

(4) 


Here,  and  throughout  this  document,  ji  is  the  gravitational  constant  GM  of  the  central  body. 
The  true  longitude  L  is  then  given  by: 


sin  L 
cos  L 


(1  —  k2b)  sin  F  +  hkb  cos  F  —  h 
1  —  h  sin  F  —  k  cos  F 
(1  —  h2b )  cos  F  +  hkb  sin  F  —  k 
1  —  h  sin  F  —  k  cos  F 


(5) 


The  third  step  is  to  compute  the  position  and  velocity  components  (. X ,  Y)  and  (X,  Y) 
of  the  satellite  in  the  equinoctial  reference  frame.  The  radial  distance  of  the  satellite  is  given 
by: 

a(l  —  h2  —  k 21 

r  =  a(l  —  h  sm  b  —  k  cos  7  )  = 


1  +  h  sin  L  +  k  cos  L 

The  position  components  are  then  given  by: 

X  =  a[(l  —  h2b )  cos  F  +  hkb  sin  F  —  k]  —  r  cos  L 
Y  =  a[(l  —  k2b )  sin  F  +  hkb  cos  F  —  h\  =  r  sin  L 

The  velocity  components  are  then  given  by: 


X  = 

Y  = 


na~ 


na~ 


[hkb  cos  F  —  (1  —  h2b )  sin  F]  = 


[(1  —  k2b )  cos  F  —  hkb  sin  F]  = 


na[h  +  sin  L) 

~  VI  -  h2  -  P 

na(k  +  cos  L) 

VI  -  h2  -  F 


Here  dots  denote  differentiation  with  respect  to  the  time  t. 

The  final  step  is  now  to  compute  the  position  and  velocity  vectors: 


(6) 


(7) 


(8) 


r  = 


r 


Xi  +  Kg 
Xi  +  Kg 


(9) 


2.1.5  Conversion  from  Position  and  Velocity  to  Equinoctial  Elements 

The  first  step  in  converting  from  position  r  and  velocity  r  is  to  compute  the  semimajor  axis 
a,  which  is  obtained  by  inverting  the  well-known  energy  integral  for  the  two-body  problem: 


The  second  step  is  to  compute  the  basis  vectors  (f,  g,  w)  of  the  equinoctial  reference 
frame.  The  w  vector  is  obtained  by  normalizing  the  angular  momentum  vector: 


Equinoctial  elements  p  and  q  are  then  given  by: 


P  = 


1  +  Iwz 


1  +  Iwz 

Vectors  f  and  g  are  then  computed  using  elements  p  and  q  in  equations  (2.1.4-la,  lb). 

The  third  step  is  to  compute  the  eccentricity-related  quantities.  The  eccentricity  vector 
e  is  given  by: 

e=-  r  rxfrxr)  (4) 


Equinoctial  elements  h  and  k  are  then  given  by: 

h  =  eg 
k  =  e  •  f 


The  final  step  is  to  compute  the  mean  longitude  A.  First  compute  the  position  coordinates 
of  the  satellite  in  the  equinoctial  reference  frame: 


X  =  rf 
Y  =  r  ■  g 


Then  compute  the  eccentric  longitude  F: 


sin  F  =  h  + 


cos  F  =  k  + 


(1  -  h2b)Y  -  hkbX 
a  a/1  -  h 2 

(1  -  k2b)X  -  hkbY 

a  a/1  -  h 2  ^k2 


where  b  is  defined  by  (2. 1.4-4).  The  mean  longitude  A  is  then  given  by  the  equinoctial  form 
(2. 1.4-2)  of  Kepler’s  equation. 
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2.1.6  Partial  Derivatives  of  Position  and  Velocity  with  Respect  to  the  Equinoc¬ 
tial  Elements 

Let 


A  =  na2  =  y/JIa 

B  =  Vl  -  h2  -  k2  =  7  -  1  (1) 

b 

C  =  1  +  p2  +  q2 


The  partial  derivatives  of  the  position  vector  r  with  respect  to  the  equinoctial  elements  are 
then  given  by: 

dr  r 


where 


da 

a 

dr 

dXn  dY 

dh 

=  - f  H - K 

dh  dh 8 

dr 

dXc  dY 

dk 

=  — f  -\ - g 

dk  dk 8 

dr 

2[Iq(Yi  -  Xg)  -  Xw] 

dp 

C 

dr 

2/[p(Xg  -  yf)  +  Yw] 

dq 

c 

dr 

r 

d\ 

n 

dX 

kX  aYY 

—  1 

~dh 

n{l  +  B)  AB 

dY 

kY  aXY 

~dh 

n(  1  +  B)  AB 

dX 

hX  aYX 

—  1  f] 

dk 

n(l  +  B)  AB  " 

dY 

hY  aXX 

dk 

n(  1  +  B)  AB 

(2) 


(3) 


The  partial  derivatives  of  the  velocity  vector  r 


with  respect  to  the  equinoctial  elements 
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are  given  by: 


dr 

da 

dr 

dh 

dr 

dk 

dr 

dp 

dr 

dq 

dr 

d\ 

where 

dx_ 

~dh 

dY 

dh 

dx_ 

dk 

dY 

dk 


2  a 

dX  dY 
dh  dh * 
dX  „  dY 
dk  dk* 

2[Iq(Yf  -  Xg)  -  Xw] 
C 

2I\p(Xg  -  Yf)  +  Yw] 
C 

naz’r 


aY 2  A  /  akX  Y2\ 

~AB  +  ^  \1  +  B  ~  ~B  ) 
aXY  A  /  okY  XY\ 

ZhT  +  ^  \TTb  +  “eT/ 

aXY  A  /  ahX  XY\ 
AB  ~  ^  \1+ ~B  +  “eT/ 
aX2  A  /  ahY  X2\ 
~~AB  ~  ^  \lTB  ~  ~B  ) 


(4) 


(5) 


2.1.7  Partial  Derivatives  of  Equinoctial  Elements  with  Respect  to  Position  and 
Velocity 

Let  and  represent  the  vectors  whose  components  in  the  (x,  y,  z)  system  are  the  partial 
derivatives  of  the  element  a,  with  respect  to  ( x,y,z )  and  (x,y,z),  respectively; 


da.i 

dr 

daj 

dr 


dai  dai  dai 
dx  dy  dz 
dai  dai  dai 
dx  dy  dz 


(1) 
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The  partial  derivatives  of  the  equinoctial  elements  with  respect  to  position  are  then  given 
by 


da 

2a2r 

dr 

dh 

abhBr 

1 

k(pX  —  IqY)  w 

B  dr 

dr 

r3  ^ 

AB 

~  ~Adk 

dk 

abkBr 

h(pX  —  IqY)  w 

B  dr 

1 _ 

dr 

AB 

a  dh 

dp 

CYw 

dr 

2  AB 

dq 

CXw 

dr 

2  AB 

dX 

r  (pX 

—  IqY)  w  bB 

( .  dr  dr 

h - h  k — 

dr 

A 

AB  A 

l  dh  dk 

The  partial  derivatives  of  the  equinoctial  elements  with  respect  to  velocity  are  given  by: 


da  2r 

dr  n2a 

dh  _  (2XY  -  AT)f  -  XXg  k{IqY  -  pX) w 
dr  n  AB 

dk  _  (2XY  —  XY)g  —  YYi  h(IqY  -  pX) w 

dr  fi  AB 

dp  CYw 

di  =  2  AB 

dq  ICXw 

di  =  2  AB 

d\  2r  k§-h§  (IqY-pX)  w 

dr  A+  1  +  B  +  A 


2.1.8  Poisson  Brackets 

The  Poisson  brackets  of  the  element  set  (oi, . . . ,  ae)  are  defined  by  the  equations 

,  ,  dai  daj  dai  daj 

°j  dr  dr  dr  dr 


(1) 


It  is  immediately  evident  that 


Oj)  0 

aj d'i) 


(2) 
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The  fifteen  independent  Poisson  brackets  for  the  equinoctial  element  set  (a,  h,  k,  p,  q,  A) 
are  given  by 


(a,  h )  =  0 

C h ,  A) 

(a,  k)  =  0 

(k,  p) 

(a,  p)  =  0 

(k,  q) 

(a,  q)  =  0 

(k,  A) 

2 

A)  =  ~a 

ip,  q) 

(ft,  ft)  =  -f 

ip,  A) 

kpC 

2  AB 

n  \  k(lC 

(h’“)=  2  AB 

(q,  x) 

hB 

A(l  +  B) 
hpC 
2  AB 
hqC 
2  AB 
kB 

A(l  +  B) 

-°L, 

4  AB 
pC 
2  AB 
qC_ 

2  AB 


(3) 


2.1.9  Direction  Cosines  (a,/?, 7) 

The  conservative  perturbations  are  more  conveniently  described  by  the  direction  cosines 
(a,/?, 7)  of  the  symmetry  axis  rather  than  the  equinoctial  elements  p  and  q.  For  central- 
body  gravitational  spherical  harmonics,  let  Zb  be  the  unit  vector  from  the  center  of  mass 
to  the  geographic  north  pole  of  the  central-body.  For  third-body  point  mass  effects  and 
shadowless  solar  radiation  pressure,  let  be  the  unit  vector  from  the  center  of  mass  to 
the  third-body.  The  direction  cosines  of  z b  with  respect  to  the  equinoctial  reference  frame 
(f,g,w)  are  then  given  by 

a  =  zb  ■  f 

P  =  zs‘g  (!) 

7  =  zb  •  w 

The  quantities  (a,/3, 7)  are  not  independent  but  related  by  the  equation 

a2  +  (32  +  72  =  1  (2) 

Note  that  (a,  f3, 7)  are  functions  of  p  and  q,  since  the  unit  vectors  (f,  g,  w)  are  functions 
of  p  and  q  through  equations  (2. 1.4-1).  Note  also  that  (a,/3, 7)  are  functions  of  t,  since  z b 
is  a  varying  function  of  time.  If  the  vector  z b  along  the  geographic  axis  of  the  central-body 
is  parallel  at  epoch  to  the  2:- axis  of  Figures  1  and  2,  then  the  direction  cosines  of  z b  are  at 
epoch 
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a 


2  Ip 


1  +  p2  +  q2 


r  1  +  p2  +  q2 

(1  —  p2  —  q2)I 

1  +  p2  +  q2 

The  partial  derivatives  of  (a,/3, 7)  with  respect  to  p  and  q  are: 


da 

2(V  +  7) 

dp 

C 

da 

2  Ip/3 

dq 

C 

d(3 

21  qa 

dp 

C 

d(3 

21  (pa  —  7) 

dq 

C 

2  a 

dp 

~C 

<9y 

21/3 

dq 

C 

(3) 


(4) 


2.2  Variation-of-Parameters  (VOP)  Equations  of  Motion 

The  Cartesian  equations  of  motion  for  an  artificial  satellite  in  an  inertial  coordinate  system 
are  [Battin,  1987]: 

i:  =  -j^  +  q  +  V77  (1) 

Here  r  is  the  position  vector  from  the  center  of  mass  of  the  central  body  to  the  satellite, 
r  =  is  the  acceleration  vector,  /j  =  GM  is  the  gravitational  constant  of  the  central 

body,  q  is  the  acceleration  due  to  nonconservative  perturbing  forces  (atmospheric  drag, 
solar  radiation  pressure),  and  77  is  a  potential-like  function  called  the  disturbing  function 
from  which  one  can  derive  the  acceleration  due  to  conservative  perturbing  forces  (central- 
body  spherical  harmonics,  third-body  point-mass).  If  m  and  n  are  the  mass  and  potential 
energy,  respectively,  of  the  satellite,  then  the  disturbing  function  77  is  given  by: 


m 


(2) 


In  order  to  apply  the  generalized  method  of  averaging,  it  is  necessary  to  convert  the 
equations  of  motion  into  a  form  giving  the  rates  of  change  of  the  satellite  orbit  elements 
as  a  function  of  the  orbit  elements  themselves.  The  equations  of  motion  resulting  from 
this  conversion  are  called  the  Variation-of-Parameters  (VOP)  equations  of  motion.  The 
derivation  of  these  equations  is  discussed  in  some  detail  in  [Cefola,  Long,  and  Holloway, 
1974]  and  [McClain,  1977]. 
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In  this  section  we  let  (oq, . .  .,a6)  =  (a,  h,  k ,  p ,  q,  A)  denote  the  osculating  equinoctial 
elements  (even  though  they  do  not  have  hats).  Then  the  VOP  equations  of  motion  turn  out 
to  be 


.  _  da.i  JU,  ,<977 

Cli  H-  ’  Q.  / 


(3) 


Here  r  =  ^  is  the  velocity  vector,  n  =  y  ^  is  the  Kepler  mean  motion,  and  <5*6  is  the 
Kronecker  delta.  The  partial  derivatives  dai/dr  are  given  by  equations  (2. 1.7-3),  and  the 
Poisson  brackets  (a*,  aj)  are  given  by  equations  (2. 1.8-2, 3). 

The  VOP  equations  of  motion  (3)  include  three  contributions  to  the  orbit  element  rates 
of  change.  The  two-body  part  is: 

a,i  =  n5ie  (4) 

The  Gaussian  or  nonconservative  part  is: 


da{ 

a*  =  -57-  •  q 
or 


(5) 


The  Lagrangian  or  conservative  part  is: 


hj 


6 


j= 1 


(6) 


In  the  remainder  of  this  document  it  will  be  convenient  to  discuss  these  contributions  sepa¬ 
rately,  but  they  must  be  added  together  to  obtain  the  total  orbit  element  rates  of  change. 

The  Lagrangian  part  of  the  VOP  equations  of  motion  contains  the  partial  derivatives  of 
the  disturbing  function  77  with  respect  to  p  and  q.  The  perturbations  which  contribute  to  7 1 
are  not  conveniently  described  in  terms  of  p  and  g,  however.  For  these  functions,  it  is  better 
to  write  7 1  as  a  function  of  (a,  h,  k,  A)  and  the  direction  cosines  ( a ,  /?,  7)  of  the  symmetry 
axis  of  the  perturbation.  The  partial  derivatives  of  the  disturbing  function  77  with  respect 
to  p  and  q  can  then  be  obtained  by  applying  the  Chain  Rule: 

<977  <977  da  <977  <9/7  <977  <9y 

dp  da  dp  +  d/3  dp  +  <9q  dp 

<977  <977  da  dll  d(3  <977  <9y 

dq  da  dq  ^  d/3  dq  <9y  dq 


The  partial  derivatives  of  ( a ,  /3,  7)  with  respect  to  p  and  q  are  given  by  (2. 1.9-4).  To  simplify 
the  notation,  let  us  again  use  the  auxiliary  quantities  A,  B,  C  defined  by  (2. 1.6-1).  Also,  let 
us  define  the  cross- derivative  operator 


_  dll  dll 

n^-a~dp  13  fa 


(8) 


Note  that  77 ,ap  =  —H,pQ.  Then  the  partial  derivatives  of  77  with  respect  to  p  and  q  turn  out 
to  be 


dll  2  /  _  _  \ 

Qp  ~  yklia'y  AlqUiap  J 

<977  2 TV  ^  \ 

q  —  (J  +p77>a/3  J 


(9) 
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With  this  notation,  equations  (6)  for  the  Lagrangian  part  of  the  VOP  equations  of  motion 
become: 


a 

h 

k 

P 

Q 

A 


2a  dB 

~Ad A 

B  dB  k 

+ 


A  dk  AB 
BdB 


IqR,,Py 


h 


Adh  AB 


~h  .  D  (  p'R'ia 7  lQ.'R'i/3'y  "b 


hB  dll 
A{l  +  B)~d\ 
kB  on 


A(l  +  B)  d\\ 


C 

2  AB 
C 

2  AB 
2a  dU 


Q  ^Tlfhk  'klia.fi 
B 

+ 


dU' 

'~d\/ 

dn' 

~d\ 


Tllft'y 

—  IB 

dn 


A  da  A(  1  +  B ) 


'  dB  , 

1  dh  +  k  dk 


+ 


AB 


(10) 


2.3  Equations  of  Averaging 

The  Generalized  Method  of  Averaging  may  be  used  to  divide  the  VOP  equations  of  motion 
(2.2-3)  into  a  short-periodic  part  which  can  be  integrated  analytically  and  a  slowly-varying 
part  which  can  be  integrated  numerically  using  time  steps  several  orders  of  magnitude  longer 
than  the  time  steps  appropriate  for  integrating  the  untransformed  equations  of  motion.  The 
Generalized  Method  of  Averaging  and  other  perturbation  techniques  are  discussed  in  [Nayfeh, 
1973].  Only  a  summary  of  the  application  of  this  procedure  to  (2.2-3)  will  be  given  here. 
More  details  can  be  found  in  [Cefola,  Long,  and  Holloway,  1974],  [McClain,  1977],  [McClain, 
Long,  and  Early,  1978]  and  [Green,  1979]. 

To  apply  the  Generalized  Method  of  Averaging  we  first  assume  that  the  osculating  orbit 
elements  a*  are  related  to  a  set  of  mean  elements  a*  by  a  near-identity  transformation: 

OO 

—  ai  +  ^2  eJrji  {a,  h,  k,p,  q,  A,  t)  (1) 

3= 1 

Here  again,  the  indexed  variables  (ai,...,a6)  refer  to  the  equinoctial  orbit  elements 
(a,h,k,p,q,  A)  and  hats  distinguish  the  osculating  elements  from  the  mean  elements.  The 
quantity  represents  a  small  short-periodic  variation  of  order  j  in  element  i.  The  quantity 
e  is  called  the  “small  parameter”  and  plays  the  role  of  a  variational  parameter  in  deriving  the 
Equations  of  Averaging.  (Note  that  the  superscript  j  is  used  in  the  symbol  e-7  to  designate 
a  power  and  in  the  symbol  r]j  to  designate  an  index.) 
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The  short-periodic  variations  are  assumed  to  contain  all  of  the  high-frequency  components 
in  the  osculating  elements  al}  so  that  the  mean  elements  a.;  vary  slowly  with  time.  This 
requirement  can  be  expressed  by  the  following  two  sets  of  inequalities: 


1 

n 


1 

n 

1 

n 

dX 

dt 


da 


dt 
da j 


dt 


—  n 


•C  a 


•Cl  for  i  —  2,  3, 4,  5 


<  1 


where  n  is  the  Kepler  mean  motion,  and 


Ak+1 

dk+1a 

dtk+l 

Ak+i 

dk+1a.i 

dtk+l 

C  a 


C  1  for  i  =  2,3, 4,  5,  6 


(2) 


(3) 


where  k  is  the  order  and  A  is  the  step  size  of  the  numerical  integrator.  Inequalities  (2)  ensure 
that  second-order  effects  will  be  small,  while  inequalities  (3)  ensure  that  the  integrator  errors 
will  be  small. 

Using  the  variational  parameter  e,  we  can  write  the  osculating  Cartesian  equations  of 
motion  (2.2-1)  as 

S  =  +<(<™  (4) 

As  e  increases  from  0  to  1  the  resulting  motion  varies  smoothly  from  two-body  motion  to 
the  actual  motion.  The  osculating  VOP  equations  of  motion  (2.2-3)  then  become 


ddi 

dt 


n(d)Sie  +  e 


ddi  JA ,  „  „  .  dTZ 


(5) 


which  can  be  written  in  the  form 


ddi 

dt 


n(a)5ie  +  tFj{af  h,  k,p,  q,  \,t) 


(6) 


The  terms  ei7)  give  the  osculating  element  rates  of  change  due  to  the  perturbing  forces  as 
functions  of  the  osculating  elements. 

We  assume  the  following  form  for  the  mean  VOP  equations  of  motion: 


OO 

(a)Si e  +  J2eJAi(a,h,k,P,q,t) 


3=  i 


(7) 
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The  terms  e^Aj  give  the  mean  element  rates  of  change  due  to  the  perturbing  forces  as 
functions  of  the  mean  elements.  For  most  perturbations,  the  A\  are  independent  of  the  fast 
variable  A  as  indicated  in  (7).  For  central-body  resonant  tessera!  harmonics,  the  A{  are  also 
slowly- varying  functions  of  j A  —  m6  (see  section  3.3). 

The  osculating  rate  functions  F,  and  the  mean  rate  functions  A{  may  be  explicit  functions 
of  the  time  t  because  the  perturbing  forces  may  change  with  time  when  the  satellite  position 
and  velocity  are  held  constant,  e.  g.,  due  to  motion  of  the  Moon. 

We  now  expand  the  osculating  rate  functions  in  a  variational  power  series  about  the  mean 
elements: 


Here 


OO 


Fi(a,h,k,p,q,X,t )  =  Fi(a,h,k,p,q,X,t )  +  E  eJ//(a,  h,  k,p,  q,  A,  t) 

i= i 


/! 

fi 


E 


3= 1 


d  Ft  l 
da3  * 


6  f)  F  1  6  6 

j=l  3  z  j= 1  k= 1 


gu  1  1 


(8) 


(9) 

(10) 


Similarly,  we  expand  the  osculating  Kepler  mean  motion  about  the  mean  semimajor  axis: 

OO 

n(a )  =  n(a)  +  ^  e-'W(a)  (11) 

3=1 


Here 


N1 

N2 

N 3 


M  +  EM1! 

2  a  8  a2 


n(a) 


3  rff  15  rjlrjf 

2  a  4  a2 


35  Ml! 

16  a3 


n(a) 


(12) 

(13) 

(14) 


Having  all  the  necessary  expansions,  we  can  now  derive  the  Equations  of  Averaging. 
First,  differentiate  (1)  with  respect  to  t  and  use  (7)  to  obtain  one  expression  for  the  osculating 
element  rates.  Next,  expand  the  functions  on  the  right  side  of  (6)  using  (8)  through  (14) 
to  obtain  another  expression  for  the  osculating  rates.  Then  equate  the  two  expansions  and 
require  that  they  be  equal  for  all  values  of  e  between  0  and  1.  Since  the  powers  of  e  are 
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linearly  independent,  the  coefficients  of  must  be  equal.  Taking  j  =  1,2,3,...  yields  the 
Equations  of  Averaging  of  order  1,2,3,...  respectively: 


Al  +  Mn(0)  + 


Aj  + 


d\ 

£>nj 

d\ 

dr]3 


n(a)  + 


dt 

drjf_ 

dt 


4  +  ty«(a)  + 94 


ax 


dt. 


(15) 

(16) 

(17) 


In  the  Equations  of  Averaging  shown  above,  the  osculating  rate  functions  Fl}  the  mean 
rate  functions  Aj,  and  the  short-periodic  variations  rjj  contain  effects  due  to  many  pertur¬ 
bations.  In  order  to  obtain  practical  expressions  for  Aj  and  rjj,  it  is  convenient  to  partition 
the  Equations  of  Averaging  to  separate  the  effects  of  different  perturbations. 

The  first  step  in  partitioning  the  Equations  of  Averaging  is  to  express  the  osculating  rate 
functions  Ft  in  terms  of  the  contributions  Fia  of  the  separate  perturbations: 

tFi  =  J2  u*Fia  (18) 

a 

The  sum  is  taken  over  all  perturbations  of  interest.  The  parameters  va  are  variational 
parameters,  one  for  each  perturbation.  Each  ua  can  vary  from  0,  at  which  the  perturbation 
is  turned  off,  to  1,  at  which  the  perturbation  has  its  actual  strength.  The  partitioned 
Equations  of  Averaging  are  required  to  be  valid  for  all  values  of  the  ua. 

Substituting  equation  (18)  into  the  first-order  Equations  of  Averaging  (15)  leads  to  the 
following  expressions  for  Aj  and  rjj: 


eA! 

^  v  ^OL-^-lOL 

a 

(19) 

er1i 

a 

(20) 

Substituting  equations  (18)-(20)  into  the  second-order  Equations  of  Averaging  (16)  leads  to 
the  following  expressions  for  A?  and  77?: 

e2A2 

^  )  ^0  yo  Ajafj 
a/3 

(21) 

e\2 

'X  )  ^'a^'fdVict/3 

a/3 

(22) 

Substituting  equations  (18)-(22)  into  the  third-order  Equations  of  Averaging  (17)  leads  to 
the  following  expressions  for  A 3  and  rjf: 

e3A3  =  VaVpV'yAiafr  (23) 

a/3"i 
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(24) 


=  Y1  uavPvlV**[h 
a/3"/ 

Similar  expressions  exist  at  higher  orders.  (Remember  that  the  first  index  in  E)a,  7fQ,  r]ia , 
etc.  refers  to  the  orbit  element.) 

If  we  substitute  (18)-(20)  into  (15)  and  then  equate  the  coefficients  of  each  ua,  we  obtain 
the  partitioned  form  of  the  first-order  Equations  of  Averaging.  Similar  procedures  with  (16), 
(17),  . . .  lead  to  higher-order  equations.  The  partitioned  Equations  of  Averaging  of  arbitrary 
order  can  be  written  in  the  concise  form 

A  +  ^ra(a)  '  (n  ~  — n(a)Sie  (25) 

Explicit  expressions  for  the  Gt  functions  up  to  order  3  are: 

Gia  =  Fiot(a,  h,  k,p,  q,  A,  t)  (26) 


Giap  =  E 

3= 1 


dFf 

dan 


-%?  + 


15  Viarii(i 
8  a2 


n(a)Si6  - 


E 


t=i 


drjia 

da3 


Ap 


(27) 


G 


ia(3 7 


/ 15 

VT  a 2 


+ 


dr)ig 

da.j 


35 1710/71/31717^  t 
16  ) n(a)Sx 


(28) 


Comparing  the  partitioned  Equations  of  Averaging  (25)~(28)  with  the  full  Equations  of 
Averaging  (15)-(17),  we  see  that  the  first-order  equations  are  identical.  The  second-order  and 
higher-order  partitioned  equations  include  auto-coupling  equations  (e.g.,  a  —  /3  —  ...  —  1) 
which  are  identical  to  the  full  equations,  but  in  addition  include  cross-coupling  equations 
(e.g.,  a  =  1,(3  =  2). 

The  partitioned  Equations  of  Averaging  (25)-(28)  give  the  fundamental  relations  which 
can  be  used  to  derive  expressions  for  the  mean  element  rates  Aia,  Aiap,  Aiat 37  and  the  short- 
periodic  variations  r)ia,  ryag ,  Then  the  total  mean  element  rates  Aj,  A2,  Af  and  short- 

periodic  variations  rj\,  rff ,  rjf  can  be  obtained  from  the  decomposition  relations  (19)-(24). 


2.4  Averaged  Equations  of  Motion 

The  Equations  of  Averaging  (2.3-25)  can  be  solved  for  the  mean  element  rates 
Aia ,  Atag ,  AjQ;37  . . .  by  applying  an  averaging  operator  <  •••  >,  to  be  defined  in  this  sec¬ 
tion,  to  both  sides  of  each  equation.  The  resulting  expressions  for  the  Aia,  Aiap,  A^^, . . . 
can  then  be  added  together  as  shown  in  equations  (2.3-19,  21,  23)  and  substituted  into 
equations  (2.3-7)  with  e  =  ua  =  1  to  form  the  mean,  or  averaged,  equations  of  motion: 

— —  =  n(a)5i6  +  E  Aa  +  E  Aap  +  E  Aap 7  +  •  •  •  (1) 

0 i  a.,/3  q:,/3,7 
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These  equations  can  then  be  integrated  with  a  numerical  integrator  to  obtain  values  for  the 
mean  elements  a*  at  a  given  time. 

The  averaging  operator  is  required  to  be  linear;  that  is,  if  p  and  cr  are  any  two  real 
numbers  and  /  and  g  are  any  two  real  piecewise  continuous  functions  of  the  mean  elements, 
then: 

<  pf  +  eg  >=  P  <  f  >  +v  <  g  >  (2) 

The  averaging  operator  is  also  required  to  be  idempotent;  that  is,  if  /  is  any  real  piecewise 
continuous  function  of  the  mean  elements,  then: 

«  /  »=<  /  >  (3) 

To  make  the  averaging  transformation  useful,  we  also  require  the  averaging  operator  to  have 
the  property  that  solving  the  Equations  of  Averaging  with  it  yields  slowly-varying  mean 
element  rates  and  small  short-periodic  variations. 

In  order  to  be  able  to  divide  the  Equations  of  Averaging  into  separate  equations  for  the 
mean  element  rates  and  the  short-periodic  variations,  we  impose  the  following  conditions: 


<  Vi  >=  0  (4) 

It  is  not  immediately  obvious  from  the  Equations  of  Averaging  (2.3-25)  that  the  short- 
periodic  variations  can  be  required  to  average  to  zero  (equations  (4)).  Let  us  first  observe 
from  (2.3-26,  27,  28)  that,  at  any  order,  Gi  are  predetermined  functions  of  the  osculating 
rate  functions  Fia  and  the  solutions  of  the  lower-order  Equations  of  Averaging.  Therefore  Gi 
are  fixed,  while  we  can  vary  and  rp  in  any  manner  which  satisfies  (2.3-25).  Let  us  assume 
that  the  short-periodic  variations  rg  do  not  average  to  zero,  and  write 


<  r)i  >=ki 


(5) 


Let  us  then  define 


. ,  .  dkt  dh  3n,  r 

Ai  =  Ai  +  n— —  +  —  +  -—kiSi6 
a  A  at  2  a 


Vi  =  Vi 


kj 


(6) 

(7) 


Solving  (6)-(7)  for  At  and  rg  and  substituting  the  resulting  expressions  into  (2.3-25)  yields 


(8) 


We  see  that  A'  and  rj[  are  solutions  to  the  Equations  of  Averaging,  and  we  can  thus  choose 
A[  and  rft  to  be  the  preferred  solutions.  Averaging  (7)  and  applying  (2,  3,  5)  yields 


<  Vi  >=  0 


(9) 


We  can  therefore  require  the  short-periodic  variations  to  average  to  zero  (equations  (4)). 
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If  the  osculating  rate  functions  Fia  for  a  given  perturbation  are  small,  27r-periodic  in  the 
satellite  mean  longitude  A,  and  slowly-varying  in  time  when  the  satellite  orbit  elements  are 
held  fixed,  then  the  single- averaging  operator  has  the  required  properties: 

<  /  >=  j  f(ai  h,  k,p,  q,  A,  t)dX  (10) 

Most  of  the  perturbations  commonly  acting  on  a  satellite  can  be  single-averaged: 

1.  Central-body  gravitational  zonal  harmonics. 

2.  Third-body  gravitational  point  mass  effects. 

3.  Atmospheric  drag. 

4.  Solar  radiation  pressure. 

Some  perturbations  are  quickly-varying  when  expressed  as  a  function  of  time  but  slowly- 
varying  when  expressed  as  a  function  of  both  time  and  a  perturbing-body  phase  angle  9 
which  varies  linearly  with  time.  If  the  osculating  rate  functions  Fia  for  such  a  perturbation 
are  small,  27r-periodic  in  both  A  and  9,  and  slowly-varying  in  time  when  6  and  the  satellite 
orbit  elements  are  held  fixed,  then  the  double-averaging  operator  has  the  required  properties: 

</>  =  ^2  J  J_  f(a,h,k,p,q,\,9,t)d\d0  (11) 

+—  cos(jX  —  m9)  f  f  f(a,h,k,p,q,X',9',t)cos(jX'^m9')dX'd9' 

2?r  0',m)eB  L  “/-7r 

+  sm(jX  —  m9)  f  f  f(a,h,k,p,q,X,,9,,t)sm(jX,  —  jm9,)dX,d9' 

J  —  IT  J  —IT 

(Equation  (11)  may  be  written  in  alternate  forms,  one  of  which  is  used  in  (3.3-2).)  Here  B 
is  the  set  of  all  ordered  pairs  (j,  m)  with  the  following  properties: 

m  >  1  (12) 

27T 

\jX  —  m9 1  <  —  (13) 

r 

where  r  is  the  minimum  period  desired  for  perturbations  included  in  the  averaged  equations 
of  motion.  The  same  minimum  period  should  be  used  for  all  double-averaged  perturbations. 
Inequality  (13)  is  called  the  resonance  condition  and  denotes  the  satellite  frequencies  j  and 
perturbing-body  frequencies  m  which  are  in  resonance  with  each  other.  It  often  happens 
that  the  satellite  has  no  resonances  with  the  perturbing  body,  in  which  case  set  B  is  empty. 
The  minimum  period  r  should  obey  the  following  inequalities: 

t  >  3  rA 
r  >  3  re 
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(14) 

(15) 


where  T\  and  tq  are  the  periods  of  the  satellite  and  the  perturbing-body  phase  angle  9 
respectively.  In  addition,  r  should  usually  be  at  least  8  times  as  long  as  the  step  size  used 
in  the  numerical  integrator,  and  much  longer  if  accurate  integration  of  perturbations  with 
this  period  is  desired.  It  is  dangerous  to  make  r  too  long,  however.  Components  of  the 
double-averaged  perturbation  which  have  periods  equal  to  r  or  shorter  will  be  excluded 
from  the  averaged  equations  of  motion  and  treated  as  short-periodic  variations.  If  r  is  too 
long  and  deep  resonance  exists,  some  of  these  short-periodic  variations  may  be  large  enough 
to  cause  large  second-order  coupling  effects,  making  the  averaging  expansions  (2.3-1)  and 
(2.3-7)  diverge.  For  this  reason,  r  should  be  less  than  100  times  the  integrator  step  size. 
Making  r  small  enough  to  ensure  convergence  of  the  averaging  expansions  takes  priority 
over  inequality  (15),  which  must  be  dropped  if  the  perturbing-body  phase  angle  9  varies  too 
slowly.  If  9  is  the  rotation  angle  of  the  Earth,  this  will  not  be  necessary,  since  the  Earth 
rotates  quickly.  If  9  is  the  rotation  angle  of  Mercury  or  Venus,  it  may  be  necessary  to  drop 
condition  (15),  depending  on  how  large  the  m-daily  (j  =  0)  short-periodic  variations  clue  to 
the  gravitational  tessera!  harmonics  are.  The  Moon  is  a  borderline  case. 

The  following  perturbations  can  be  double-averaged: 

1.  Central-body  gravitational  sectoral  and  tessera!  harmonics.  For  this  perturbation,  9  is 
the  rotation  angle  of  the  central  body.  (For  a  more  precise  definition  of  9,  see  Section 
2.7.1)  The  double-averaged  central-body  gravitational  spherical  harmonic  model  in 
SST  is  fully-developed  and  will  be  discussed  in  Section  3.3. 

2.  Third-body  point-mass  effects,  if  the  third  body  orbits  the  central  body.  For  this 
perturbation,  9  is  the  equinoctial  mean  longitude  of  the  third  body.  The  current  double- 
averaged  third-body  model  assumes  that  the  third  body  orbits  the  central  body  in  a 
slowly-varying  Keplerian  ellipse.  Methods  for  predicting  the  effects  of  short-periodic 
variations  in  the  third-body  orbit  on  the  satellite  orbit  have  yet  to  be  developed.  For 
Earth  satellites,  the  short-periodic  variations  in  the  orbit  of  the  Moon  are  substantial. 
If  they  are  included  in  the  Lunar  ephemeris  used  with  the  double-averaged  third-body 
model,  the  integrator  step  size  will  be  driven  down  to  values  appropriate  for  the  single- 
averaged  third-body  model,  thus  destroying  the  usefulness  of  the  double-averaging 
expansion.  The  step  size  reduction  is  avoided  by  using  a  smoothed  ephemeris  for  the 
Moon,  but  this  creates  an  error  in  the  computed  satellite  mean  element  rates  clue  to 
the  Lunar  perturbation,  and  the  size  of  this  error  is  not  known.  Because  of  these 
limitations,  double-averaged  third-body  perturbation  models  will  not  be  discussed  in 
detail  in  this  document.  For  a  complete  description  of  the  current  model,  see  [Collins, 
1981], 

There  are  some  perturbations  for  which  no  averaging  operator  with  the  required  proper¬ 
ties  can  be  found.  These  perturbations  are  called  non-averageable  and  include: 

1.  Atmospheric  drag,  with  an  asymmetric  spacecraft  and  fast,  non-periodic  attitude  vari¬ 
ation. 

2.  Solar  radiation  pressure,  with  an  asymmetric  spacecraft  and  fast,  non-periodic  attitude 
variation. 
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3.  Continuous  thrust,  with  fast,  non-periodic  changes  in  direction. 

4.  Impulsive  thrust. 

These  perturbations  are  typical  of  directed  flight,  which  in  general  is  not  required  to  be  either 
slowly-varying  or  27r-periodic  in  A  or  any  other  phase  angle.  Semianalytic  Satellite  Theory 
cannot  predict  the  effects  of  these  perturbations,  with  the  exception  of  impulsive  thrust  (see 
below),  and  should  not  be  used  unless  they  are  small  enough  to  ignore. 

Some  types  of  directed  flight  are  averageable,  and  these  include  many  scenarios  of  prac¬ 
tical  interest: 

1.  A  drag-perturbed  satellite  whose  solar  panels  always  point  directly  at  the  Sun.  The 
attitude  of  the  satellite  relative  to  the  atmosphere  will  be  27r-periodic  in  A  and  will 
vary  slowly  in  time  as  the  Earth  revolves  about  the  Sun. 

2.  A  spacecraft  with  a  solar  sail  which  is  feathered  when  approaching  the  Sun  and  per¬ 
pendicular  to  the  Sun  line  when  receding.  The  attitude  of  the  sail  will  be  27r-periodic 
in  A  and  will  vary  slowly  in  time  as  the  Earth  revolves  about  the  Sun. 

3.  A  spacecraft  with  a  constant-thrust  ion  engine  whose  thrust  is  always  parallel  to  the 
orbit.  The  direction  of  the  resulting  acceleration  will  be  27T-periodic  in  A  and  the 
magnitude  will  vary  slowly  in  time  as  reaction  mass  is  depleted  and  the  spacecraft  gets 

lighter. 

Of  course,  the  perturbations  remain  averageable  only  as  long  as  the  orbit  remains  elliptical. 
Parabolic  and  hyperbolic  orbits  are  beyond  the  scope  of  this  document. 

Impulsive  thrust  is  not  averageable,  but  its  effects  can  be  predicted  using  the  following 
procedure: 

1.  Integrate  the  averaged  equations  of  motion  (1)  up  to  the  time  of  the  impulsive  maneu¬ 
ver. 

2.  Compute  the  short-periodic  variations  as  functions  of  the  mean  elements  using  equa¬ 
tions  (2.3-20,  22,  24)  and  the  equations  in  Section  2.5. 

3.  Add  the  short-periodic  variations  to  the  mean  elements  to  get  the  osculating  elements 
(equations  (2.3-1)). 

4.  Convert  the  osculating  equinoctial  orbit  elements  to  position  and  velocity  using  the 
equations  in  Section  2.1.4. 

5.  Add  the  velocity  change  Av  caused  by  the  impulsive  maneuver  to  the  satellite  velocity. 

6.  Convert  the  satellite  position  and  velocity  to  osculating  equinoctial  orbit  elements  using 
the  equations  in  Section  2.1.5. 

7.  Invert  equations  (2.3-1)  to  convert  the  osculating  elements  to  mean  elements  (see  Sec¬ 
tion  6). 
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This  procedure  is  always  valid,  but  if  impulsive  thrusts  occur  more  often  than  once  per  orbit, 
it  may  be  too  expensive  to  make  the  use  of  averaging  worthwhile.  For  additional  methods 
for  modeling  impulsive  maneuvers,  as  well  as  continuous  thrust,  see  McClain  [1982], 

Some  perturbations  are  usually  averageable,  but  cannot  be  averaged  in  certain  circum¬ 
stances  because  of  large  second-order  effects.  These  include: 

1.  Third-body  point-mass  effects,  for  a  satellite  whose  orbit  comes  too  close  to  the  orbit 
of  the  third  body.  This  perturbation  will  be  non-averageable  even  if  resonance  locking 
ensures  that  the  satellite  will  always  remain  far  from  the  third  body. 

2.  Third-body  point-mass  effects,  for  a  satellite  whose  orbit  comes  too  close  to  the  bound¬ 
ary  of  the  central  body’s  gravitational  Sphere  of  Influence. 

3.  Atmospheric  drag  effects,  during  the  terminal  stage  of  reentry. 


Semianalytical  Satellite  Theory  should  not  be  used  under  these  circumstances. 

It  is  worth  considering  at  this  point  whether  the  averaging  operators  (10)  and  (11)  ac¬ 
tually  have  the  required  properties  (2)-(3).  It  is  immediately  clear  from  equations  (10)  and 
(11)  that  both  operators  are  linear  (equation  (2)).  It  is  also  clear  from  equation  (10)  that  the 
single-averaging  operator  is  idempotent  (equation  (3)).  To  show  that  the  double-averaging 
operator  is  idempotent,  we  double-average  both  sides  of  equation  (11)  to  obtain 


«/» 


<  f  >  d\  dO 


2i T2 


Y  [cos(jA  —  m6 ) 


+  sin(jA 


<  /  >  cos(jA' 


f  <  f  >  sm(jX' 

J  —7T 


-  md')d\'dd' 

-  m9')d\' d6'] 


1 

47T2 


fdXdO 


/  —7 r  J  —  7T 


+— -  Y,  [cos(j\  —  m9)  f  f  cos2 (j X' —  m9')dX'd9'  f  f  f  cos(jX''  —  m9")dX!'d,9" 

+  sin(jA  —  m9)  f  f  sin2(j'A/  —  m9')dX!d9'  f  j  f  sm(jX"  —  m9")dX"d9"] 

J  —IX  j  —IT  J  —TV  J  —IT 

=  </> 

(16) 

Note  that,  if  /  is  independent  of  9 ,  the  double-averaging  operator  (11)  reduces  to  the 
single-averaging  operator  (10).  Thus  we  can  without  inconsistencies  apply  the  single-averaging 
operator  (10)  to  perturbations  that  do  not  depend  upon  9,  and  the  double-averaging  operator 
(11)  to  perturbations  that  do  depend  upon  9. 

If  the  perturbations  can  be  averaged  with  the  operators  (10)  or  (11),  then  the  Equations 
of  Averaging  (2.3-25)  can  be  solved  for  the  mean  element  rates  by  averaging  both  sides  of 
each  equation.  Note  that  the  averages  are  always  taken  with  (a,  h,  k,p,  q,t)  held  constant. 
Using  (4),  we  thereby  obtain  for  the  mean  element  rates  A*  at  any  order 


Ai  — <  Gi  > 


(17) 
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Substituting  (2.3-26,  27,  28)  into  (17),  we  can  obtain  explicit  expressions  for  the  mean 
element  rates  at  each  order. 

The  explicit  equations  for  the  first-order  mean  element  rates  have  the  same  form  for  both 
single-averaged  and  double-averaged  perturbations: 

Aia  =<  Fia(a,h,k,p,q,X,t)  >  (18) 


The  averaging  operator  used  is  single  in  the  first  case  and  double  in  the  second. 

The  explicit  equations  for  the  second-order  mean  element  rates  involve  coupling  terms 
between  two  perturbations,  which  may  be  two  different  perturbations  (cross-coupling)  or 
copies  of  the  same  perturbation  (auto-coupling).  If  either  perturbation  is  single-averaged, 
the  equations  take  the  form 


15  Tl 

H  3  2  (VlaVip)  ^*6 

8  cr 


(19) 


If  both  perturbations  are  single-averaged,  a  single- averaging  operator  is  used.  If  either 
perturbation  is  double-averaged,  a  double-averaging  operator  is  used.  If  both  perturbations 
are  double-averaged,  the  equations  take  the  form 


A 


ia/3 


.  15  n  .  /  ^Via 


A?/3 


(20) 


If  both  perturbations  use  the  same  perturbing-body  phase  angle  9,  then  a  double-averaging 
operator  is  used.  However,  if  the  perturbations  use  different  perturbing-body  phase  angles 
9  and  6',  then  a  triple-averaging  operator  is  used. 

Derivation  of  the  explicit  equations  for  the  third-order  mean  element  rates  will  be  left  as 
an  exercise  for  the  reader. 


2.5  Short-Periodic  Variations 

Once  the  mean  element  rates  Aia,  Aiag,  Aiapy, . . .  are  known,  the  Equations  of  Averaging 
(2.3-25)  can  be  solved  for  the  short-periodic  variations  Via ,  rjiag ,  Viapy<> ...  by  using  Fourier 
series  expansions  for  the  functions  Gia,Giap,Giapy, . . ..  The  resulting  expressions  for  the 
Via,  Via/3,  Viap-y,  ■  ■  ■  can  then  be  added  together  as  shown  in  equations  (2.3-20,  22,  24)  with 
e  =  ua  =  1  and  substituted  into  equations  (2.3-1)  with  e  =  1  to  give  the  osculating  elements: 

Oj  0>i  T  )  )  Via  T  )  )  Viafl  T  )  )  Via/3 7  T  •  •  •  (1) 

&  a/3  a(3~/ 

Of  course,  using  Fourier  series  expansions  for  the  functions  G*  assumes  that  the  osculating 
rate  functions  Fia  are  27r-periodic  in  the  phase  angles  of  the  expansions.  Most  perturbations 
can  be  expressed  with  more  than  one  kind  of  Fourier  series  expansion.  In  the  following 
subsections  we  shall  give  formulas  for  the  short-periodic  variations  corresponding  to  several 
possible  expansions. 
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Some  of  the  general  expressions  presented  in  this  section  are  new,  so  enough  steps  in 
the  derivations  will  be  presented  to  enable  the  reader  to  rederive  them.  Our  derivations 
are  based  on  the  work  of  [Cefola  and  McClain,  1978],  [Green,  1979],  [McClain  and  Slutsky, 
1980],  [Slutsky,  1980],  [Slutsky  and  McClain,  1981],  and  [McClain,  1982], 


2.5.1  General  rjt  Expansions  in  A 

Recall  that  the  Equations  of  Averaging  of  arbitrary  order  can  be  written  as 


A 


drji  drji  3  n 

+  XT  =  Gi  ~  o -0*6?7 1 
cm  at  2  a 


where  Gi  are  predetermined  functions  (equations  (2.3-26,  27,  28))  and 

A  =<  Gi  > 


(1) 

(2) 


At  order  m,  the  Equations  of  Averaging  contain  m  perturbations,  which  may  all  be  different 
or  may  include  multiple  copies  of  the  same  perturbation. 

If  all  of  the  perturbations  in  the  Equations  of  Averaging  (1)  are  single-averaged,  then  the 
equations  can  be  solved  for  the  short-periodic  variations  rjt  by  integrating  over  the  satellite 
mean  longitude  A.  It  is  convenient  to  first  define  the  short-periodic  kernels 


6  =  -  [ (Gi  -  Ai)d\ 
n  J 

The  constant  of  integration  is  specified  by  requiring 


<  &  >=  0 


(3) 

(4) 


In  the  absence  of  explicit  time-dependence,  the  full  short-periodic  variations  r)i  are  then  given 
by: 


(5) 


The  conditions  (2.4-4)  and  (4)  require 


/ 


=  0 


(6) 


In  the  presence  of  explicit  time-dependence,  the  full  short-periodic  variations  77*  are  given 
by: 


Vi 


K  (— lA 

e,  +  E(  j 

k=  1 


nK 


'2  Ji6 


d% 

dtk 


K 


dX 


£idA  +  +  1) 


(-1)* 


dkG 


k= 1 


n 


k  Jk+ 1  dtk 


d\k+i 


(7) 


Here  K  is  the  order  of  the  highest  order  partial  derivatives  desired  in  the  expansions.  The 
conditions  (2.4-4)  and  (4)  require 


0  for  all  k  >  1 


(8) 
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Here  the  following  somewhat  unusual  notation  is  used  for  the  operator  which  is  the  kth 
indefinite  integral  (the  inverse  of  the  kth  derivative  operator): 


J  f(X)d\t  =  j  ■■■  J  /(A)  dA'-VA. 


(9) 


Alternately,  knowing  the  short-periodic  variations  77°  in  the  absence  of  explicit  time- 
dependence,  we  can  compute  recursively  the  short-periodic  variations  rjk  including  the  kth 
order  time  derivatives: 

ri°  =  (1-  [  (idX 

“  >  3k  r  +„»  (10) 

1 wd>f+^-^L  h 


t  =  t-i  .  ("It  [3**1° 

li  li  '  k 

nr  Jk 


-dxk+1 


4+1  dtk 

Explicitly,  we  suppose  that  the  functions  G+  can  be  written  as  a  Fourier  series  in  the 
mean  longitude  A: 

Gi(a,  h,  k,p,  q,  A,  t)  =  C°(a,  h,  k,p,  q,  t) 


(11) 


+  X  [Ci  (a>  4  k,  P,  q,  t)  cos  jA  +  SI  (a,  h,  k ,  p,  q,  t )  sin  j  A 
3= 1 


where 


C?  =  — 


G,  dX 


2n  J- 71 
1  />7r 

C{  =  —  GiCosjXdX 

71  J-ir 

Sj  =  —  I  Gt  sin  j  X  dX 

7T  J  —  7T 

Using  (2.4-10)  and  (2),  we  can  obtain  the  mean  element  rates  Ap 

A,  =  C? 


(12) 


(13) 


We  can  then  obtain  the  short-periodic  variations  p,  by  performing  the  integrals  (7,  8).  The 
final  result  of  these  calculations  is: 

OO 

Vi  =  X  i°i  cos  j X  +  Si  sin  J  M  ( 14) 

j= i 


where 


ci 

si 


3  1  dig  j 

2aJC' 


+ 


U 


n) 


dC[  32  S^dSi 
dt  2  a  j  dt 


1 

1 

r+y 

3  3  4  d2C{' 

1 

' d3C{ 

3  4  4,  4341 

1 

CO 

_  dt2 

2  a  j  dt 2 

Cm)4 

_  dt3 

2  a  j  dt3 

3  1  4>  j 


+ 


(j 


n 


;-~Si 

■a  j 

2C{  3  3  5i6d2S{ 

tt2  2a  j  dt 2 


/  L 

'dSj  _  3  2Si6dC[ 


dt  2a  j  dt 

1  rg<g  _  344,  <93Cj 

(jn)4  <9t3  2a  j  dt3 


(15) 


2.5.2  General  Tjt  Expansions  in  F 

In  this  section  we  suppose  that  the  functions  Gi  are  expanded  as  a  finite  Fourier  series  in 
the  eccentric  longitude  F: 

Gi(a,h,k,p,q,F,t)  =  Cf(a,h,k,p,q,t )  (1) 

j 

+  X  lCi  (a’ /z’  k’  t)  cosJF  +  sl  (a,  h,  k,  p,  q,  t)  sin  jF 

3= 1 

where 

c“  =  hSlG'iF 

C\  =  —  r  Gi  cos  jF  dF  (2) 

7 r  J —ir 

S{  —  —  [  Gi  sin  j  F  dF 

7 r  J —n 

Again  we  may  use  (2.4-10)  and  (2.5. 1-2,  3,  4,  7,  8)  to  obtain  the  short  periodic  variations 
rji.  To  do  this,  it  will  be  necessary  to  convert  the  integrals  over  A  into  Fourier  series  expansions 
over  F.  We  begin  by  supposing  /(A)  has  a  finite  Fourier  series  expansion  in  F  with  known 
coefficients  and  it  averages  to  zero: 


/(A)  =  c°  +  J2(cJ  cos  JF  +  & sin  :iF) 

3= 1 


<  /(A)  >=  0 


Using 


J  f{X)d\  —  j  f(X)—dF 

the  following  consequence  of  the  equinoctial  form  (2. 1.4-2)  of  Kepler’s  Equation 

dX  r 

— —  =  -  =  1  —  h  sin  F  —  k  cos  F 
oF  a 

and  the  following  well-known  trigonometric  identities 

cos(j  —  k)F  +  cos  (j  +  k)F 


cos  jF  cos  kF  = 
cos  jF  sin  kF  = 
sin  j  F  cos  kF  = 
sin  jF  sin  kF  = 


sin  (j  +  k)F  —  sin(j  —  k)F 
2 

sm(j  +  k)F  +  sin(j  —  k)F 
2 

cos(j  —  k)F  —  cos  (j  +  k)F 


(3) 

(4) 

(5) 

(6) 


(7) 


29 


we  can  convert  the  right  side  of  (5)  into  a  Fourier  series  expansion  in  F.  Note  in  particular 
that  the  condition  (4)  implies  that  the  constant  term  C°  in  (3)  must  be  related  to  the  Fourier 
coefficients  C1  and  51  by 

c°  =  +  is1  (8) 

Higher-order  integrals  can  be  computed  using  recursion  formulas  obtained  from  the  equa¬ 


tion 


[  f(X)dXm+1  =  [  I  f(X)dXmdX 

Jm-\- 1  J  Jm 


(9) 


We  summarize  the  conversion  for  the  general  multiple  integral  with  the  following  notation: 


J+m 


[  f(X)dXm  =  U°m(C f ,  S<)  +  £  [Ui(C^  &)  cos  jF  +  V' (<*  &)  sin  jF]  (10) 

3= 1 


<  f  f(X)d,Xm  >=  0 

J  m 


in) 


Here  the  arguments  denote  that  the  functions  UL  and  Vi  depend  upon  the  coeffi¬ 

cients  and  appearing  in  (3)  through  the  relations: 
for  1  <  j  <  J\ 

ui  =  v 


for  m  >  0 


Vo  =  & 


Ul  = 


- u 1  +  -V1 

2  m  1  2  171 


TJ 1 

u  m+1 


— U 1  -  (l  -  — \  V1  -  -U2  +  -V2 

2  m  i  2  J  rn  2  m  2  m 


for  2  <  j  <  J  +  m  +  1: 


=  1 


Ul  -  —Vi 


-u2  -  -V2 

2  m  2  m 


(12) 


JP 

um+ 1 


J 


h  k  h  h 

-vj  +  -tH'-'  +  -Vj~l  -  -Uj+1  +  -W+1 

r  m.  1  ^  m  1  r  m  ^  m  1  r  m 


C 


for  j  >  J  +  m  +  1: 


=  -  ( Uj  -  - U j~l  +  - Vj~l  -  -Uj+1  -  -Vj+1 

■  \  m  2  m  2  m  2  m  2  m 

v  \  / 


UL  =  0 


vz  =  o 
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With  this  notation,  we  can  write  the  mean  element  rates  Ai  and  the  short-periodic 
variations  rji  corresponding  to  (1): 


where 


Cl 

si 


A  =  c°  -  - C 1  -  -S1 

1  1  Tl  Tl 

J+K+2 

Vi  =  C°+  £  \ci  COSJF  + Si  sm  jF 
j= i 


C»  =  k-cl  +  h-sl 


n 


u;  (cf.sf)  +  y. 


-i)' 


uL 


I 


dk6  dkSl 


) 


'2  Ji6 


A  nk+l  k+l  \  3tk  ’  8tk  ) 

&  (d.tf)  +  £(* + (A-  if 


fc=i 


»„(««) 


n 


VI  (cf,5{)  +  Yik  + 


(13) 

(14) 


(15) 


2.5.3  General  r/?:  Expansions  in  L 

In  this  section  we  suppose  that  the  functions  Gj  are  expanded  as  a  finite  modified  Fourier 
series  in  the  true  longitude  L: 


Gi(a,  h,  k,p,  g,  L,  t) 


M 

Ci(a,  h,  k,  p,  q,t)+J2  Am(«,  h,  k,  p,  q,  t )  (L  -  X)m 

m=  1 
J 

+£[b  (a,  h ,  k,p,  q,  t )  cos  jL  +  5/  (a,  h ,  k,p,  q,  t )  sin  jL 
i=i 


(1) 


Here  the  quantities  (L  — A)m  are  written  separately,  rather  than  by  replacing  them  with  their 
Fourier  series  expansions  recorded  below. 

The  equation  of  the  center  may  be  calculated  from  the  Fourier  series  expansion 

oo  2 

L«-  A  =  V]  -(<7j  cos  jL  —  pj  sin  jL)  (2) 

j=i  J 


Here 


1  />7r 

Pj  =  (cos  jL)  —  —  /  cos  jLdX 

ZlT  J  —  7T 

1  />7r 

cr,  =  (sin  jL)  —  —  sinjLdA 

27T  J-tt 


(3) 
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The  auxiliary  quantities  p3  and  a3  can  be  evaluated  using  the  equations  * 

Pj  =  (1  +jB)(-b)jCj(k,h) 

=  (1  +  3B){-b)3Sj{k,h) 

where  b  and  B  are  given  by  (2. 1.4-4)  and  (2.1.6-lb),  and 

Cj(k,  h)  +  iSj(k,  h)  =  (k  +  ih)j 
are  obtained  from  the  recursion  formulas 


Cj+1{k,  h )  =  kCj(k,  h )  -  hSj(k,  h)  ,  C0  =  1 
Sj+±(k,  h )  =  hCj(k,  h)  +  kSj{k ,  h)  ,  Sq  =  0 
The  quantities  (L  —  \)m  may  be  calculated  from  the  expansion 

OO 

(L  -  \)m  =  n°m  +  ]T«.  cos  jL  +  sin  jL) 

3= 1 

Here  the  first-order  coefficients  k{  and  ipl  are  obtained  immediately  from  (2): 

re?  =  0 

4  = 

W  = 

J 


(4) 

(5) 

(6) 

(7) 

(8) 


Higher-order  coefficients  K,Jm  and  E  are  obtained  using  the  following  recursion  formulas, 
obtained  by  multiplying  the  series  on  the  right  sides  of  (2)  and  (7): 

OO  ^ 

«m+l  =  E  ~  PfcV’m) 

fc=l  K 


j  2a j  0  ~  1 

J  fc=i  7  + 

oo  i  i-i  i 

+  E  T(^<+fe  -  +  (1  -  <5,l)  E  —^-(°3-kKkm  + 

k= i  K  fe=i  4  K 

V'm+l  =  E-VT^i+^m+Pl+^m) 

3  k= 1  3  ' 

oo  i  j-1  1 

+  E  7 ~S<Tki’mk  +  +  (1  —  <5jl)  E  “  r(Cri-fc'0m  ~  Pj-k^m) 

k=  1  K  fc=l  4  K 


"Let  Z  =  exp(?X)  and  use  the  Residue  Theorem  to  integrate 

(1  —  h2  —  /c2)3/2  expfijX) 


p3  +  i(Tj  — 


27T 


/_w  (1  +  h  sin  L  +  k  cos  L): 


-gX 


(9) 
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Note  that  if  m  is  odd,  then  (L  —  A)m  is  antisymmetric  about  perigee.  Therefore  we  have 

1 


K2k+1 


(10) 


(L  —  A  )2i+1dL  =  0 

27 r  J —7T 

i.  e.  ,  only  even  powers  of  (L  —  A)  have  nonzero  constant  terms. 

Again  we  may  use  (2.4-10)  and  (2.5. 1-2, 3, 4, 7, 8)  to  obtain  the  short-periodic  variations 
rjt .  To  do  this,  it  will  be  necessary  to  convert  the  integrals  over  A  into  modified  Fourier  series 
expansions  in  L.  We  begin  by  supposing  /(A)  has  a  modified  Fourier  series  expansion  in  L 
with  known  coefficients  and  it  averages  to  zero: 

M  j 

/(A)  =  C°  +  E  v'n(L  -  x)m  +  E(ci  cos  JL  +  sj  sin iL) 

3= 1 


m=  1 


<  /(A)  >=  0 


Using 


r  c  \ 

j  f(X)d\  =  J  f(X)—dL  , 

the  following  consequences  of  (2. 1.4-2, 5, 6)  and  (2) 

dX  1  /r\2  (1  —  h2  —  A;2)3/2 


(11) 

(12) 

(13) 


dL  VI  -  h 2  -  k2  \a 
the  integral 


(l+hsinL  +  kcosLY  =^Uf>i^jL  +  ajsinjL),  (14) 


( T  —  A')m+1 

L  -  X)mdX  =  -1 - 1 - +  /  (L-  X )mdL  , 

m+  1 


(15) 


the  identities  (2. 5. 2-7)  and  expansions  (7),  we  can  convert  the  right  side  of  (13)  into  a  Fourier 
series  expansion  in  L.  Note  in  particular  that  equations  (3)  and  (12)  imply  that  the  constant 
term  C°  in  (11)  must  be  related  to  the  Fourier  coefficients  Cl,S‘  and  T>1  by 


M  J 

c°  =  -E®”«S.-E(cV4  +  5' 

m=  1  j= 1 


<7i 


(16) 


Note  also  that  we  need  the  following  formulas  for  the  product  of  two  Fourier  series: 


J  OO  1  oo 

E(£J  cos  j  L  +  S3  sin  jL)  E(Pfc  cos  kL  +  Ok  sin  kL)  =  -  N{  l((j)(C3pJ+S3 

3= 1  '  ?  '  ' 


k= 1 


3= 1 


(1  -  E)  E  (V-'Pt  -  sj~e°e)  +  Ji/_1(j)  eV+%  +  sj+%) 

f=max(j-  J,l)  £=1 

J 


■  +  SC<Tj+e) 


e=  i 


cos  jL 


(i  -  E)  E  E(-ci+V<  +  sj+epe) 

1=1 
J 


E  (CV3+,  -  5'ft+() 


£=1 


sin  j  L  | 


(17) 
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Here  Z*(j)  is  the  inclusion  operator  defined  by 


1  if  r  <  j  <  s 
0  otherwise 


(18) 


Higher-order  integrals  can  be  computed  using  recursion  formulas  obtained  from  equation 
(2. 5. 2-9).  We  summarize  the  conversion  for  the  general  multiple  integral  with  the  following 
notation: 


f(X)d\k 


M+k 

V  WT(C<,S<,TV)(L-Ar 

m=  1 
oo 

+  y  \ui^,S(,V<)cosjL  +  Vi^,S(,V<)SmjL 

3= 1 


(19) 


Qj(\)d\k}=0  (20) 

Here  the  arguments  (tb  ,  S ‘q  V '»)  denote  that  the  functions  Uk,  V£,  and  W™  depend  upon  the 
coefficients  C\S^,  and  T>m  appearing  in  (11)  through  the  relations 
for  j  >  1: 


U(  =  — 


M 


3- 1 


vi  = 


1 

j 


+  E  +  (1  -  ^i)  E  ^~kak  +  Si-kpk) 

m=  1  fc=max(j-  J,l) 

E(-^+vfc  +  sj+kPk)  +  E(cVJ+fc  -  skPj+k) 

k= 1  k= 1 

M  j- 1 

i/(i)c'  +  y  r<  +  (i  -  a,,)  y  (c>-yt  -  s‘-+k) 

m=l  fc=max(j- J,l) 

•t-J  J 

E(ci+V*  +  ^+vfc)  +  E(cVfc  +  sv,+fe) 


fc=i 


fc=i 


H-V  =  -C° 

/jym— 1 


for  k  >  1: 


m 


M+k 


ui  =  -  y  w-rc  -  y<dft- + he 

ra=l  j= 1 


d+1  = 


M+fc-1 


1-1 


h+  y  + (i  -  <^0  y (ufat + vfp,) 


m=  1 


1=1 


+ E(-^  ^  +  tfe*  +  vrPi  -  nV*) 

€=1 


(21) 
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v;+ 1  = 


M+k-1 


3- 1 


E  + (i E(^rV( - ^r'^) 


m=l 


e= i 


+  E(^+V  +  £4W  +  ^J+E  +  nS+*) 

1=1 


Wlk+ 1  =  -c/£ 


W&1 


m 


With  this  notation,  we  can  write  the  mean  element  rates  A j  and  the  short-periodic 
variations  rjt  corresponding  to  (1): 


Ji  aj ) 


where 
C?  = 


M  J 

A  =  C\  +  E  V? +  iZiCjpj  +  S{' 

m= 1  j=l 

M+if+2  oo 

I?i  =  c“+  X  Drti-Ar  +  ^C/cosjL  +  SfsinjL) 

ra=l  j=l 


M+K+2  oo 

x  Am<  -  E(dft  +  s U,) 

m=  1  j= 1 


(22) 

(23) 


Ci  =  w<  5<  ,  f-  wiv 

+ A,  nt+i  t,*+i  ^  Qtk  •  g,»  •  9(1-  J 


n 


A' 


-h;, 


20'“  tE2  s?  ’  ®f> + e<* + yy,  yy) 


vAc;,SiS®9  +  E 

*  =  /5*cf  ^  a‘^ 

k= 

K 


n 


^  ^ j  _  '-'j  _  ‘-'i 

rj  nfc+1  fc+1  V  <9ffc  ’  <9tfc  ’  dtk  ) 


(24) 


2a5 

D™  =  Zf+1 


r,  i  n(~1)tr.i  jggj  gd  ggh 


2o.«l-Ki^,5j,®f)  +  2rl*+  DW^dlt'ir-arj 


m 


-irriQ.^-X?)  +  (^2, 

n  "  1  v  ’  nk+1  k+l  \  <9tfc  <9ffc  <9ffc  / 


_3_a 

'2a'5"’ 


K 

-E2? 

*:=! 


J  i  A _ EL_ 


{m)(k  +  1) 


Mt™  /avf  (ASf  a^f\ 


n 


fc+1  ^fc+2  ^  ’  ftfc  ’  dtk  ) 


In  the  absence  of  explicit  time-dependence,  equations  (23)-(24)  can  be  simplified.  The 
short-periodic  variations  become: 


M+ 2 


m -C?+£  DT(L  ~  XT  +  E(C?  cos  ji  +  s?  sin  ji) 

ra=l  j  =  l 


(25) 
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where 


c?  =  -  £  °t<  -  Zicipi  +  si”,) 

771=1  j= 1 

c{  =  -u;(cls},v<)-^-sieuUci,s<,vi) 

n  Zan 

S’  =  ivy(Cf,  5f,  Cf)  -  5^-iieV?'^,  Sf,  P<) 

n  2an  (26) 

A1  =  --c,°  +  Wfec/ftef,  sf,  ®{) 

n  2  an 

Di  =  -—vl  ~  t —<*«£? 

1  2n  1  4 an  1 

1  s 

oh  i\  ®i§D\ 
mn  2am(m  —  l)n 

2.5.4  General  rjt  Expansions  in  A,  0 

If  one  or  more  perturbations  are  double-averaged,  then  the  functions  Gt  can  be  written  as  a 
double  Fourier  series  in  the  mean  longitude  A  and  the  perturbing-body  phase  angle  Q: 

Gi(a,  h,  k,p,  q,  X,  9,  t)  =  y ^\Cjm(a,  h,  k,p,  q,  t )  cos(jX-mO)  +Sjm(a,  h,  k,p,  q,  t )  sin(jA  —  m6)\ 


where 


c00  = 


47T2  J — 7r 


GidXdO 


Cjm  =  — —  f  /  Gt  cos(jA  —  m6)d\d6 

Z7T  J  —  7T  j  —TV 


gjm  _  x 


Gi  sin(jA  —  mQ)d\dQ 


27T2  7-tt  7-tt  ' 

Using  (2.4-11)  and  (2. 5. 1-2),  we  can  obtain  the  mean  element  rates  Ap 

A  =  ^  [C/m  cos(jA  —  mQ)  +  S[m  sin(jA  —  m0)]  (3) 

C 3,m)eB 

We  can  then  obtain  the  short-periodic  variations  77*  by  integrating  the  Equations  of  Averaging 
(2.5. 1-1)  in  a  manner  similar  to  the  single-averaged  case  in  Section  2.5.1.  Assuming  that  the 
coefficients  Cjm,  S(rn  and  the  rotation  rate  6  do  not  explicitly  depend  upon  time,  we  obtain 

rji  =  ]T)  [Cjm  cos(jA  —  m6 )  +  Sfm  sin(jA  —  m6)]  (4) 


where 


(J3  m  _ 


gjm  = 


jn  —  mQ 
1 

- r  I 

jn  —  mQ 


2  a  jn  —  mQ 

Qjm  ^  3  n  hjg 


2  a  jn  —  mQ 
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2.5.5  First-Order  T)ia  for  Conservative  Perturbations 


For  conservative  perturbations,  it  may  be  advantageous  to  use  an  alternate  solution  for  the 
first-order  short-periodic  variations  rjia  which  avoids  having  to  obtain  the  osculating  rate 
functions  Fia .  If  a  perturbation  a  is  conservative,  then  the  osculating  rate  functions  Fia  can 
be  expressed  as  (2.2-6): 

ID 


3  da  j 


Fia  'y  ' ,  iflii  Q'i 
3= 1 

where  (aq . . .  a§)  are  the  equinoctial  elements  (a,  h,  k,p ,  g,  A),  the  quantities  (a*,  dj )  are  the 
Poisson  brackets  given  by  (2. 1.8-2, 3),  and  1Z  is  the  osculating  disturbing  function. 

If  the  perturbation  a  is  single-averaged,  we  can  define  the  mean  disturbing  function  U : 

1  rn 

U=<7Z>=—  /  7Z(a,  h,  k,p,  q,  X,t)dX  (2) 

27T  J-ir 

Averaging  both  sides  of  (1),  we  obtain  the  first-order  mean  element  rates  due  to  this  pertur¬ 
bation: 

5  dU 

■Aia  y  ) ((Z j ,  Qj )  —  (3) 

j= 1  aa3 

The  first-order  short-periodic  variations  caused  by  this  perturbation  can  be  obtained  from 
a  potential-like  function  S  called  the  short-periodic  generating  function: 


S  =  J(U-  U)dX 


(4) 


<  S  >  =  0 


(5) 


Using  (2.3-26),  (2.4-10),  and  (l)-(4),  we  can  obtain  the  first-order  short-periodic  kernels  6 
by  performing  the  integrals  (2. 5. 1-2,  3)  (with  the  subscript  i  replaced  by  ia ): 


,  ,dS_ 

—  Z^va*’  a3)  an 

U  -  =  l  UUj 

From  (5)  it  is  clear  that: 

(6a)  =  0 

Combining  (2. 1.8-2, 3)  and  (6)-(7),  we  obtain: 

2  ds 


6a 


n2a  dX 
2 


/  6«dA  =  —  S 
J  nza 


(6) 

(7) 

(8) 

(9) 


Substituting  (6)  and  (9)  into  (2.5. 1-5),  we  obtain  the  following  expression  for  the  first-order 
short-periodic  variations  r)io  in  the  absence  of  explicit  time-dependence: 


V  ioc 

n 


4-3  3dS  3  c  c 

/  j ( O'i ;  Q j )  ^  T  o^iS^ 

1=1 


da.j  na2 


(10) 
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In  the  presence  of  explicit  time-dependence,  the  first-order  short-periodic  variations  rjlo  are 
given  by: 


i  J.f,  ,as 
n>°  „ 


K 

E- 

k= 1 


-i  Y 


TV 


5  8  r  3k  S' 

E(a*>  ai)"»T"  I.  TuiT^k  +  (a*> 


3  e 

d~  ^  0*6 
na- 


■=1  7  aaj  Jk  dtk 

S  +  J2(k  +  l)^^-  f  ^d\k 

n  A  dtk 


r  dW 

L-  i  dffc 


dA 


fc-i 


in) 


2.5.6  Second-Order  r/,Qi3  for  Two  Perturbations  Expanded  in  A 

In  this  section  we  suppose  that  the  osculating  rate  functions  iyQ  and  the  first-order  short- 
periodic  variations  7/ja  can  be  written  as  Fourier  series  in  the  mean  longitude  A: 

Fia(a,h,k,p,q,  A)  =  C°a(a,h,k,p,q,t) 

OO 

+  E[<E(a>  h,  k,p,  q,  t )  cos  jA  +  S/a(a,  h,  k,p,  q,  t )  sin  jA]  (1) 
j=i 


Via  =  Y(CiacosJX  +  SL  sinjA) 
i=i 

From  (2.3-27),  the  second-order  functions  G*a/t  are 


(2) 


dF. 


15n , 


^-y  kj±IOL  .  x  ^ 6  r-  V— V  ^  •  \lOL  a 

Giaf)  =  §  WE/r/9  +  8^67?larh/3  “  §  dE4-'3 


(3) 


Substituting  (l)-(2)  into  (3),  and  using  (2.5.3-17)  with  J  =  oo,  we  can  write  Gja/3  as  a 
Fourier  series  in  A: 

OO 

GW  =  C°a/3  +  Y,  (cL/3  cos  Jx  +  sLp  sin  3- x)  (4) 

3= 1 

where 


F°  =  -  V 

'-'ia/3  q  /  j 


5  /  OyJ  O  03  \  1  c.,n 

E|  ^ia  qJ  \  |  _  ,;/7j  ct  i  A__A.  trd  rd  J-  Q-?  Q-?  \ 


2  E  Lr=l  V  da 


8  a2 


CW  =  E 


r=l 


dCf 


-G;'  -C 


dCi  1  -  8,jX  W  /dc£rfc 


0 


da.  r/3  r/3  da. 


+ 


E 


'^20'  s~ik 


CZn  ~ 


dSl 


-k 

ioi  nk 


2  WiV  dar  rP  dar  r/3 


1  00  f  f)r]+k  r)Ck  8^+k  N 

2  E  r/3  dW  r/3  "dW  fc/3  dW  fc/3  , 


1  _  A  i_1  r  1  E..V, 

-iCSeSL  +  — ^ E(i- *0(Si  C?„  +  C%)  +  g^MC^c*, - s^%) 


fc=l 
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-|  OG 

E  [0  +  k)(si:kc%  -  c£ksk„)  +  k(slc$k  -  eisif) 


k= 1 


i  s:  fr<3+krik  ,  nk  r<j+k  ,  o3+kok  ,  ok  o3+k\ 

'  8a2  °i6 ^°la  "T  Dla  Dl/3  +  d1oD1/3  J 


5L„  =  V 


r=l 


rr?  _  /i0  w^icx. 

~dZ  r0~  r0~dZ 


dsL  .  i-<V^Vdt±"E,  dstk 


+ 


E 


ia  0&;  i  ^  *  o  ^jk 


-<%  + 


2  ^  V  <9ar  r/3  Qar  r0 


_lg/  dC£k  dd 

2  k= i  ' 


Cfc  ,  dCiaoi+k  .  dSl  rk  dSLro+k\ 

dar  br0  +  <9ar  ^  +  <9ar  <9ar  ^  J 


1  -  5,1  ^ 


fc=i 


15n 


+WL  +  —7^  E  (J  -  k)(Sl~kSk6l3  -  Ci~kCk6(3)  +  —Si6  (C£%  +  Sjrcy 


(5) 


-I  CXJ 

+5  E  [-«  +  *)(4,+%  +  c;:kc *,)  +  k(s*,sl f  +  c‘ c&‘) 

z  k= 1 

i  (  rij+kak  i  ^-yfc  ci+fc  i  ctf+fc/^fc  cfc  /^b+fcN 

"t"  g  2  \  L  la  Dl/3  '  l'1o°1/3  +  Dlo  °l/3  Dla(-"i/3  ) 

Then  the  second-order  mean  element  rates  are 

A  _  /'’O 

^ija/3  —  L'ja/3 

and  the  second-order  short-periodic  variations  are 


(6) 


Via0  =  E(C'ECOS^A  +  Sia0sln:jX) 
3=1 


(7) 


where  the  Cjag  and  SJiag  are  given  by  (2.5.1-15)  (with  the  subscript  i  replaced  by  ia/3)  in 
terms  of  Cjag  and  Sjag.  The  formulas  in  this  section  were  published  in  [Danielson,  March 
1993], 


2.5.7  Second-Order  pia^ 3  for  Two  Perturbations  Expanded  in  L 

In  this  section  we  suppose  that  the  osculating  rate  functions  Fia  and  the  first-order  short- 
periodic  variations  r)ia  can  be  written  as  finite  modified  Fourier  series  in  the  true 
longitude  L : 

Fia(a,h,k,p,q,L,t)  =  C^(a,h,k,p,q,t ) 

Ja  (1) 
+  EI  CL(a,h,k,p,q,t)  cos  jL  +  Sja(a,  h,  k,p ,  q,  t)  sin  j L] 

3=1 
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(2) 


Via(a,h,k,p,q,L,t)  =  C^(a,h,k,p,q,t)  +  53  DZ(a,  h,k,p,q,t)(L  -  X)m 

771=1 

Ka 

+  53  [Cia(a,  h,  k,  p,  q,  t )  cos kL  +  Ska(a,  h,  k,  p,  q,  t )  sin  kL] 


k= 1 


The  second-order  functions  Gja/3  are  again  given  by  (2. 5. 6-3).  Differentiating  (1)— (2) ,  we 
can  obtain  the  needed  partials 


9Frn  =  CZ  +  53  {C&  cos  jL  +  SZ  sin  jL) 

3  = 1 

Mar  Kar 

=  C£  +  E  D.T(i  -  A)"‘  +  E  (c£  cos  kL  +  St  sin  kL) 

m=  1  k= 1 


<9ar 

drpc 

dar 


(3) 

(4) 


The  product  of  two  Fourier  series  can  be  converted  into  a  single  Fourier  series  with  the 
formula 

J  K 

53  cos  jL  +  S-7  sin  jL)  ^2(Ck  cos  kL  +  Sk  sin  kL) 

j= 1  k= 1 

i  J+A' 

=  -  X]  {xrn(J’/°(j)(CiC'J  +  SLS”) 

Z  3= 1 

min(jf  —  1  ,-RT)  min(J— j,K) 

+  [J/+A'(j)  ^  (Ci_fcCfc-5J'“fc5,fc)+J/-1(j)  X]  (CJ'+fcCfc  +  5J'+fc5fc) 

/c=max(j  — J,l)  fc=l 

min(if— j,  J) 

Tjf-1  (j)  53  (CfcCJ+fc  +  5fe^'+fe)]  cos  jL 

fc=l 

min(j— 1,K)  min  (J—j,K) 

+[x i+K{j)  53  (cj~ksk  +  sj-kck)  +  i^-\j)  53  (-<x+fesfe  +  sj+kck) 

/c=max(j  — J,  1)  /c=l 

min^— j,  J) 


+l?~1(j)  53  (CfcLj+fc  -  5fcCJ'+fc)]  sin  jL} 


(5) 


fc=i 


With  the  use  of  (2. 5. 3-7),  we  then  obtain 


M“r+M/3  oo 

G.«s  =  C+  E  ®Ss(i-Ar  +  E(CEcosjL  +  5EsinjL)  (6) 

m=l  j= 1 


40 


where 


C°  =  V 

'~'iaf5  /  j 


r=  1 


Jar  Ml3 

czc% +we  (C4, + si:rjD?0 
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1 5n 
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^ ia(3 
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z  fc=l 
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+^r-(i)  e  (c£ci+t + 5£si+t) 

z  fc=i 

MP 

-ir^ArfiCt+T,  D^rU)&m 
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+5(1  -  E  (<&-*"•<&  - 

fc=max(j- Jar,l) 

J“r-j  J“r 
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k= 1 
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15  n 


(j)cici, + 1 row. 
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^sa{ii“‘v)cisi u+irojc^ 
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r=l 

1 5n 


m— 1 


+  8a2^6 


xrH^Q^+xrH^rQ+xr+M/3H  e 

j=i 


1/3 


Then  the  second-order  mean  element  rates  are 

Mar+M^  oo 

A-iaf}  =  Cia/3  +  E  C/3Km  +  E(^ial8Pi  +  ^L/3aj)  (§) 

m=l  j=l 

and  the  second-order  short-periodic  variations  are 

M“r+Ar^+A'+2  oo 

>fa»  =  <?»,„  +  E  ^(i-ir  +  EfC^cosiL  +  S^sinii)  (9) 

m=l  j=l 

where  the  C\ap  ,  SJia/3  ,  D™g  are  given  by  (2.5.3-24)  (with  the  subscript  i  replaced  by  ia/3 
and  M  replaced  by  Mar  -f-ilX3)  in  terms  of  C3ia*  and  SjQ/3.  The  formulas  in  this  section  were 
published  in  [Danielson,  August  1993]. 


2.6  Partial  Derivatives  for  State  Estimation 

Observational  data  may  be  used  to  improve  the  estimate  of  a  satellite’s  state.  Some  differ¬ 
ential  correction  algorithms  which  have  been  used  in  conjunction  with  SST  are  described  in 
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[Green,  1979],  [Taylor,  1982]  and  [Long,  Capellari,  Velez,  and  Fuchs,  1989].  It  is  our  purpose 
here  only  to  explain  how  to  obtain  the  partial  derivatives  needed  in  such  Liters. 

We  let  Ok  denote  the  value  of  the  kth  observed  quantity  computed  with  the  SST  orbital 
generator.  The  SST  state  variables  are  the  initial  mean  elements  Gq(to)  and  various  constant 
parameters  q  (the  geopotential  coefficients  Cnm  and  Snrn  in  (2.7-1),  the  drag  coefficient  Co 
in  (3.4-3),  the  solar  radiation  pressure  coefficient  Cr  in  (3.5-6),  etc.  ).  Required  for  a  batch 
filter  are  the  partial  derivatives  of  the  Ok  with  respect  to  the  state  variables  a*(fo)  and  q. 

The  actual  observations  are  commonly  of  position  and  velocity  components  in  a  local 
coordinate  frame  fixed  on  the  surface  of  the  Earth.  However,  through  transformations  these 
components  may  be  expressed  in  terms  of  the  orbital  elements,  so  we  can  regard  the  Ok  as  an 
implicit  function  of  the  osculating  elements  dj.  Application  of  the  chain  rule  then  produces 

d Ok  A  dOk  daj 

d a* (t0)  ddj  ddiito ) 

dCh  _  ®  d(hda^ 

dci  ddj  dci 

Assuming  we  can  obtain  the  partials  analytically,  our  remaining  task  is  to  calculate  the 
partials  y  and  Differentiating  the  decomposition  (1-1)  yields 


da. 


dai(t0 


k=i 


ddj  JA  / 

ULl  k=l  \ 


|  drjj  \  ddk 
ddk )  ddi{t0 ) 

(3) 

drjj  \  ddk  !  drjj 
ddk)  d^  d^ 

(4) 

The  partials  d^a(htu )  and  are  often  arranged  to  form  a  matrix  $,  referred  to  as  the 
state  transition  matrix: 


$(Mo)  = 


dd  i 

dd  i 

dd  i 

dd 

ddiito) 

dae  (to) 

dci 

dcy 

ddo 

ddo 

ddo 

dd{ 

dd^to) 

ddo  (ffi) 

dci 

dc , 

0 

0 

1  0 

...  o 

0  1 

1  0 

0 

0 

0  ••• 

0  1 

(5) 


Here  £  is  the  number  of  parameters  q.  Differentiating  (5)  with  respect  to  t  and  interchanging 
the  ordinary  and  partial  derivatives,  we  can  obtain  the  following  initial  value  problem  for 
$(Mo): 

3>(Mo)  =  F&(t,t0%  *(t0,tQ)=I  (6) 
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Here  I  is  the  identity  matrix  and 


da  1 

da  1 

0 

da  1 

da 

da  1 

da5 

dci 

dct 

dae 

dae 

0 

dciQ 

dd, 

da± 

da5 

dc\ 

dc, 

0  .  0 


(7) 


0 


0 


The  $  matrix  is  a  function  only  of  the  five  slowly  varying  mean  elements,  and  therefore 
the  numerical  integration  of  (6)  can  be  done  with  the  same  large  step  size  as  used  in  the 
integration  of  equations  (l)-(2)  for  the  mean  element  rates.  Values  of  <f>  at  observation  times 
not  coinciding  with  the  integrator  step  times  can  be  obtained  by  interpolation. 

Our  task  has  thus  been  reduced  to  the  calculation  of  the  partial  derivatives 

and  These  same  partials  are  also  needed  in  a  sequential  Kalman  filter.  For  the  two-body 
part  of  the  mean  element  rates 

ai  =  ndi6  (8) 


the  only  nonzero  partial  in  (7)  is 

dde  3  n 

da  i  2  a 


(9) 


and  thus 


$ 


1 

0 

0 

0 

0 

3  n(t  —  to) 
2  a 
0 


0 

1 


0  1 


(10) 


0 


1  0 
0  1 


Although  the  differential  correction  algorithm  for  updating  the  initial  mean  elements  may 
converge  with  only  the  two-body  partials  (10),  the  speed  of  convergence  can  be  improved 
by  including  the  dominant  perturbation  partials.  Analytical  formulas  have  been  obtained 
for  the  partial  derivatives  with  respect  to  the  mean  elements  at  of  the  J2  contribution  to 
the  mean  element  rates  a j  (equations  (3.1-12)),  the  partial  derivatives  with  respect  to  the 
geopotential  coefficients  Cnm  and  Snm  of  the  resonant  tessera!  contribution  to  the  mean 
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element  rates  a j  (equations  (3.2-9)),  and  the  partial  derivatives  with  respect  to  the  mean 
elements  ax  of  the  J2  contribution  to  the  short  periodic  variations  rj3  (from  the  expressions 
in  Section  4.1). 


2.7  Central-Body  Gravitational  Potential 

The  well-known  expression  for  the  disturbing  function  due  to  the  gravitational  held  of  the 
central  body  is  [Battin,  1987]: 

N  min (n,M)  /  n 

K{r,  (p,ip)  =  -Y\  E  (  —  )  (sin  0)  (Cnm  cos  imp  +  Snm  sin  imp)  (1) 

T  n=2  ,7^0  ^  r  7 


radial  distance  from  center  of  mass  of  central  body 

geocentric  latitude 

geographic  longitude 

central-body  gravitational  constant 

central-body  mean  equatorial  radius 

associated  Legendre  function  of  order  m  and  degree  n 

geopotential  constant  coefficients 

maximum  order  of  geopotential  held  (M  <  N) 

maximum  degree  of  geopotential  held 


In  this  section  all  elements  are  osculating  (even  though  they  do  not  have  hats). 

In  the  hrst  subsection  we  shall  outline  the  development  of  the  central-body  gravitational 
potential  into  the  form  used  in  SST.  Complete  details  are  to  be  found  in  [Cefola,  1976], 
[McClain,  1978],  [Proulx,  McClain,  Early,  and  Cefola,  1981],  and  [Proulx,  1982],  Then  in 
the  remaining  subsections  we  describe  methods  for  calculating  the  various  functions  used  in 
the  expansion. 


2.7.1  Expansion  of  the  Geopotential  in  Equinoctial  Variables 

We  start  by  writing  (2.7-1)  in  the  complex  form 

{  N  min (n,M)  /  R\n  1 

~E  E  (  — )  -Pnm (sin 0)(Cnm  -  iSnm)exp(irrnp)  >  (1) 

r  n^2  ^ 0  ^  J 

Here  i  =  \f—\  and  Re  {^}  is  the  real  part  of  z.  With  the  goal  of  expressing  (1)  in  terms 
of  the  equinoctial  elements,  we  set 

ip  =  aB  -  9  (2) 

where  6  is  the  central  body  rotation  angle  and  a#  is  the  right  ascension.  If  we  let  (xB,  Yb,  zb) 
denote  a  right-handed  orthonormal  triad  fixed  in  the  central  body,  with  xB  pointing  to  the 
prime  meridian  and  zB  to  the  geographic  north  pole,  then  9  may  be  calculated  from  [Early, 
1982]: 
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(3) 


sin  9  = 


cos  6  = 


-f  •  yB  +  /g  •  xB 

1  +  Py 

f  •  xg  +  Jg  •  yg 

1  +  Py 


(Remember  that  7  is  defined  by  (2.1.9-lc).) 

Next  the  spherical  harmonics  P„;m  (sin  0) exp  (ima:#)  are  expanded  as  a  Fourier  series  in 
the  true  longitude  L,  using  a  rotational  transformation  theorem  for  the  spherical  harmonics: 

n 

Pnm(sin  0) exp (imaB )  =  K™'STSexP(*sI/)  (5) 

s=—n 

The  VP  coefficients  are  defined  by: 


vrt  = 


(- 1 (n  +  s)!(n  —  s)! 

2n  (n-m)!(2±2)!(^)! 


if  n  —  s  is  even 


if  n  —  s  is  odd 


The  rotation  functions  S™s(a,f3, 7)  may  be  expressed  in  terms  of  the  dot  products  (a,  (3, 7) 
of  the  Zb  vector  with  the  equinoctial  reference  triad  (f.  g,w): 


gms  _  jr, 


(rl)m-*2 s(a  +  iP)Im~s(  1  +  1 7) ~Im P™+ss~m-s (7)  if  s  <  -m 

^  (nFw)!(fl  lTl)\  at  n\m—Is(-\  ,  T~,\Is  r>m—s,m+s  /  n,\  ;  (•  |„|  ^  ^ 

(-1)  2  (n  +  6.)|(n_  g)|  +  (!  +  J7)  Pn-m  (7)  if  Is!  —  m 


2~s(a  -  i0)s“/m(l  +  Ii)ImPnZTS+m(l) 


if  s  >  m 


(Remember  that  /  is  defined  by  (2. 1.2-2).)  ffere  P/w( 7)  are  Jacobi  polynomials.  (Note  that 
commas  are  used  to  separate  indices  in  a  symbol  such  as  P^Xss~m~s  in  order  to  prevent 
ambiguities.) 

Next  the  product  (r)nexp(isL)  is  expanded  in  a  Fourier  series  in  the  mean  longitude  A 
(the  sixth  equinoctial  element  a@): 


exp (isL)=  jr  YJlsexp(ij X) 


Now  the  Hansen  coefficients  X^s  are  defined  by  [Hansen,  1855] 


r\n  00 

-  j  exp  (is/)  =  £  Xf  exp(ijM) 

a  j=- 00 
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and  the  kernel  Kr-'s  of  the  Hansen  coefficient  X™s  is  defined  by 


KJs{e)  =  e~\s-j\X™(e)  (10) 

where  /  is  the  true  anomaly,  Ad  is  the  mean  anomaly,  and  in  (10)  e  is  the  orbital  eccentricity. 
Hence,  remembering  equations  (2. 1.2-1,  4),  we  can  express  the  functions  YJiS  as 

YJs(h,  k)  =  [k  +  ih  sgn  (s  -  (11) 


Here 


sgn(x) 


1  if  x  >  0 
—  1  if  x  <  0 


(12) 


The  last  step  is  simply  a  rearrangement  of  the  order  of  summation,  so  as  to  isolate  the 
total  disturbing  potential  due  to  the  phase  angle  j A  —  m9,  and  to  facilitate  the  use  of  stable 
recursion  formulas.  We  also  introduce  some  new  notation,  to  enable  the  results  to  be  written 
concisely.  First,  define  the  functions 


r„ms(7) 


Next,  put 


< 


(_  i)"*-»2s(l  +  /7)-/m 

r_1  (n  +  m)\(n  —  m)\ 

1  J  (n  +  s)!(n-s)! 


(1  +  J7)/s 


2-s(l  +  /7)/m 


if  s  <  —  m 
if  |s|  <  m 
if  s  >  m 


(13) 


Q3  i  j  fji 

m.ft  1  '  m. 


[k  +  ih  sgn(s  —  j)]ls  jI(ck  +  Uj3)m  Is  if  \s\  <  m 

[k  +  ih  sgn(s  —  j)}^s~j^[a  —  i(3  sgn(s  —  m)]^s~Im^  if  |s|  >  m 


(14) 


Then  dehne  the  Jacobi  polynomial  P™  indices  by 


£ 

v 

w 


[  n  —  m  if  |s|  <  m 

In—  |s|  if  |s|  >  m 
m  —  s\ 

m  +  s\ 


(15) 


The  disturbing  function  can  now  finally  be  written  as 


K  =  Re{^  E 


M 

E 


N 

E 


N 

E 


j=—o o  m= 0  s=—N  n=max(2,m,|s|) 


R 


jmympm  j^—n—l,Spvw 


(16) 


' 6';',, .  +  ■//;(,,]( C -  15„m)exp[l(jA  -  md )]} 


Note  that  the  functions  GJms  and  HJms  defined  by  (14)  are  of  degree  s  —  j  in  the  eccen¬ 
tricity.  The  power  s  —  j  has  been  called  a  D’Alembert  characteristic  (not  in  accordance 
with  its  original  connotation  [Brouwer,  1961]). 
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2.7.2  Calculation  of  V™  Coefficients 

The  V™  coefficients  are  defined  by  (2. 7. 1-6).  Since 


Vm  =  (-lYVm 

v  n,s  \  )  v  n,—s 


(1) 


we  can  restrict  onr  discussion  to  the  case  s  >  0  without  loss  of  generality.  Furthermore, 
since  V)™.  =  0  when  n  —  s  is  odd,  we  need  only  consider  the  case  when  n  —  s  is  even.  Also, 
note  that  the  lowest  value  of  the  degree  n  in  the  summations  (2.7.1-16)  is  greater  than  or 
equal  to  2 ,  m,  and  s | . 

Suitable  recurrence  relations  are 

(n  +  s  +  l)(n-s  +  l) 

(n  —  m  +  2){n  —  m  +  1)  n’s  (2) 

(n  -  m)V™s 

Appropriate  initialization  is  provided  by 

=  1 

T/°  =  ^±Ht/0  ^ 

S+1’S+1  (s  +  1)  s’s 

That  is,  to  calculate  the  V™  coefficients,  first  use  (3)  to  get  values  for  rn  —  0  and  n  —  s  for 
s  —  0, 1, . . .  .  Then  use  (2b)  for  m  >  0  and  still  n  =  s.  Finally,  use  (2a)  for  increasing  n 
with  any  nonnegative  m  and  s. 


Vm 

v  n+2,s 


V. 


m+1 


n,s 


2.7.3  Calculation  of  Kernels  K'J1S  of  Hansen  Coefficients 

From  the  definitions  (2. 7.1-9,  10),  the  kernels  of  the  Hansen  coefficients  are  given  by 


AT 


e  = 


g  Is  3 1  ri r 

2n  J-TT 


cos  (sf  —  jM)dM 


(1) 


The  kernels  of  the  Hansen  coefficients  are  thus  functions  of  the  orbital  eccentricity  e. 
from  (1)  that 


Note 

(2) 


so  we  can  restrict  our  discussion  to  the  case  s  >  0  without  loss  of  generality. 

For  the  special  case  j  —  0,  the  kernels  may  be  evaluated  in  a  form  algebraically  closed 
in  the  eccentricity.  This  is  because  the  Hansen  coefficients  with  j  =  0  are  related  to  the 
associated  Legendre  functions  by  [McClain,  1978] 


\r  —  n—l,S 


\rn,S 

Ao 


in-  l)\Xr‘ 


Pn-i(x)  for  n  >  1 


(n  +  s  —  1)! 

(— l)s(n  —  s  +  l)!y~n~1 
(n  +  1)! 


Pn+i(x)  for  n  >  0 


(3) 
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where 


_  p  _  k 2  ' 

and  many  recursion  formulas  are  available  for  the  associated  Legendre  functions,  which  for 
arguments  in  the  range  1  <  x  <  00  are  defined  by 

1  An+s 

Kix)  =  ^(x2  -  i)s/2^(x2  -  1)”  (5) 

For  the  special  kernels  with  the  first  superscript  negative  (needed  in  Sections  3.1  and  4.1), 
appropriate  recursion  formulas  are 


r—n—l,s 


for  n  —  s  >  0 
for  n  —  s  +  1  >  1 


(n  -  l)y2 

(n  +  s  —  1  )(n  —  s  —  1) 


"(2 n  -  3 )Kon,s  -  in  -  2)/l0“n+1’s]  for  n  >  s  +  2  >  2 


For  the  special  kernels  with  the  first  superscript  nonnegative  (needed  in  Section  3.2),  appro¬ 
priate  recursion  formulas  are 


Kns  — 


with  initializations 


(2s  —  1)  t^s_2,s-1 
-h-o 
s 

(2s  +  1)  g_i,s 
s  +  1 

2n  +  1  1)S  (n  +  s)  (n  -  s)  n_2)S 

i  1  ^0  /  |  -i  \  2  ^0 

n  +  1  n[n  +  l)x 


for  n  —  s  —  1  >  1 


for  n  —  s  >  1 


for  n  >  s  +  1  >  2 


K°0’°  =  1 

o  i  (8) 

Ko’1  =  -1 

The  general  kernels  K~n~1,s  may  be  computed  from  the  following  recurrence  relation 


[Proulx,  McClain,  Early,  and  Cefola,  1981]: 


(3  —  n)(l  —  n  +  s)(l  —  n  —  s) 


- -  {(3  —  n)(l  —  n)(3  —  2n)Kj 

—  n.  —  <?)  I*  J 


2js-\Ts-n+l,s  ,  -2/ 


-(2  -  n)[(3  -  n)(  1  -  n)  +  —  }K;n+l’s  +  j2(l  -  n)K] 


r— n+3,s 


To  initialize  this  recurrence  relation,  we  need  the  values  of  the  four  kernels  K]  n’s,  Kj  n+1,s, 
KJn+ 2,s,  and  KJn+3,s  at  n  =  max (2 ,m,s).  These  latter  kernels  are  calculated  by  infinite 
series  representations.  In  a  study  of  the  various  possibilities  [Proulx  and  McClain,  1988],  it 
was  found  that  the  expansion  of  choice  is 


Kf  =  (1  -  e2)n+i  ^ 


rns  2  a 

a-\-a,a+b^ 
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Here 


a  =  rnax(j  —  s,  0) 
b  =  max(s  —  j,  0) 


in) 


The  Y™  terms  are  called  modified  Newcomb  operators  and  may  be  computed  by  the  recur¬ 
rence  relation 


4  (P  +  a)Ypy  =  2(2  s-n)Y£tt  +  (s-n)YF&*-2(2s  +  n)YZ£11 

—  (s  +  n)Yp£_  2  +  2(2p  +  2a  +  2  +  3n)Ypl\)a_i 
Recursion  formula  (12)  is  initialized  by 


(12) 


vn,s 

I  n 


0,0 


=  1 


(13) 


and  by  treating  quantities  with  negative  subscripts  as  identically  zero.  That  is,  to  calculate 
the  Y™  coefficients  for  any  n,  use  (12)  -(13)  to  get  values  for  s  —  0, . . . ,  n  and  each  successive 
value  of  p  =  0, 1, . . .  and  a  =  0, 1, . . .  Note  that  the  Y™  terms  are  rational  constants,  and 
therefore  they  can  be  computed  once  and  stored  for  all  later  applications. 


2.7.4  Calculation  of  Jacobi  Polynomials  P 

The  Jacobi  polynomials  appear  in  the  expression  (2. 7. 1-7)  for  the  rotation  functions.  P^w(y) 
is  a  polynomial  of  degree  £  in  7,  which  from  (2.1.9-lc)  is  the  cosine  of  the  angle  between 
a  vector  from  the  geographic  south  to  north  pole  of  the  central  body  and  the  angular  mo¬ 
mentum  vector  of  the  satellite.  The  Jacobi  polynomials  P}mi  (7)  with  m  —  0  in  the  indices 
(2.7.1-15),  i.e.  i  —  n  —  s  >  0  and  v  =  w  =  s  >  0,  are  related  to  the  associated  Legendre 
functions  Pns( 7)  by 

K 1.(7)  =  2'rw-r?)1  -  T2)-/2fts(7)  (1) 

[n  +  sj! 

The  Jacobi  polynomials  can  be  computed  from  the  standard  recurrence  relation  [Szego, 
1959]: 


2£(£  +  v  +  w)(2£  +  v  +  w  -  2)P?W('Y)  = 

(2£  +  v  +  w  —  l)[(2f  +  v  +  w)(2£  +  v  +  w  —  2)7  +  v 2  —  tc2]P/JfL(  7)  (2) 

-2[£  +  v-l  ){£  +  w-  1){2£  +  v  +  w)P?™2{  7) 

This  recursion  formula  is  initialized  by 


JDvw 
r 0 


pvw 


1 

0 


(3) 
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2.7.5  Calculation  of  G3ms  and  H3ns  Polynomials 

From  the  definitions  (2.7.1-14),  the  functions  G3ms  and  H3ms  are  polynomials  in  the  equinoctial 
elements  h,  k  and  the  direction  cosines  a,  /3. 

These  polynomials  may  all  be  calculated  from  one  set  of  generic  recurrence  formulas, 
based  on  the  Cj  and  Sj  polynomials  obtained  from  (2. 5. 3-6): 


G3 

ms 


C\s-j\(k,  h)Cm_Is(a,  (3)  -  isgn(s  -  j)S\s-j\(k,  h)Sm_Is(a,  (5)  for  |s|  <  m 

C\s-j\(k,  h)C\s_Im\(a,  (5)  +  sgn(s  —  j)sgn(s  —  m)S\s-j\(k,  h)S\s_Im\(a,  (5)  for  |s|  >  m 

(1) 


H3 

ms 


IC\s-j\(k,  h)Sm-l9(a,P)  +  sgn(s  -  j)S\8-j\(k,  h)Cm_Is(a,  (3)  for  |s|  <  m 

— sgn(s  —  m)C\s-j\{k,  h)S\s-Im\(a,  (3)  +  sgn(s  —  j)S\s-j\(k,  h)C\s_Im\(a,  j3)  for  |s|  >  m 

(2) 


2.8  Third-Body  Gravitational  Potential 


The  disturbing  function  due  to  the  gravitational  field  of  a  third-body  point  mass  is  [Battin, 
1987]: 


77(r,  0,  t) 


/i3  l  R3  r  cos  0 
-R3  \  |R3  —  r|  R3 


(1) 


Here 

r  =  vector  from  the  center  of  mass  of  the  central  body  to  the  satellite 
R3(t)  =  vector  from  the  center  of  mass  of  the  central  body  to  the  third  body 
0  =  angle  between  the  vectors  r  and  R3 
fi3  =  third-body  gravitational  constant 

In  this  section  all  elements  are  osculating  (even  though  they  do  not  have  hats). 
The  quantity  |R^r|  can  be  expanded  in  the  following  series: 


Rs 


R3  -  r 


(2) 


where  Pn  is  the  Legendre  polynomial  of  degree  n.  Hence  the  third-body  disturbing  function 
(1)  can  be  written  as 

E(f)n^n(cOS0)  (3) 

K'i  n= 2 

where  N  is  the  maximum  power  of  the  parallax  factor  to  be  retained  in  the  expansion. 

The  further  development  of  the  third-body  potential  into  the  form  used  in  SST  is  similar 
to  that  of  the  central-body  potential  in  Section  2.7.  Complete  details  are  to  be  found 
in  [Cefola  and  Broucke,  1975],  [McClain,  1978],  [Cefola  and  McClain,  1978],  and  [Slutsky, 
1983], 
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2.8.1  Expansion  of  Third-Body  Potential  in  Equinoctial  Variables 

With  the  goal  of  expressing  (2.8-3)  in  terms  of  the  equinoctial  elements,  we  set 

cos  (j)  —  a  cos  L  +  (3  sin  L  (1) 


where  L  is  the  true  longitude,  and  now  ( a ,  /?,  7)  are  the  dot  products  of  the  unit  vector 
with  the  equinoctial  reference  triad  (f,  g,  w).  Due  to  the  motion  of  the  third  body, 
(a,  /3, 7)  are  slowly  varying  functions  of  the  time  t  and  are  the  source  of  weak  time-dependence 
effects. 

Next  the  expression  Pn(a cos  L  +  fj  sin  L)  is  expanded  into  a  Fourier  series,  using  an 
addition  formula  for  the  Legendre  polynomials: 

n 

Pn{a  cos  L  +  (3 sin  L)  =  ^(2  -  80s)VnsQns{^)[Cs(a,  (5)  cos  sL  +  Ss(a,/3)  sinsL]  (2) 

s=0 


Here  <5os  is  the  Kronecker  delta,  and  the  polynomials  Cs(a,/3)  and  Ss(a,/3)  are  the  same  as 
those  in  (2. 7.5-1,  2).  The  new  coefficients  introduced  in  (2)  are  defined  by: 


K, 


0 


if  n  —  s  is  even 
if  n  —  s  is  odd 


(3) 


Qns  (5/) 


dsPn(  7) 

djs 


After  substituting  (2)  into  (2.8-3),  we  can  write  the  result  in  the  complex  form 


1^3  n=2s=o 


n  /  r\n 

VnsQns[Cs(a,  P)  -  iSs(a,/3 )]  exp (isL) 


(4) 

(5) 


The  last  steps  are  simply  the  replacement  of  the  product  elsL  by  the  expansion 
(2. 7. 1-8),  and  a  rearrangement  of  the  order  of  summation.  The  disturbing  function  can  then 
be  written  as 


2.8.2  Calculation  of  Vns  Coefficients 


The  Vns  coefficients  are  defined  by  (2. 8. 1-3).  Since  Vns  =  0  when  n  —  s  is  odd,  we  need  only 
consider  the  case  when  n  —  s  is  even.  (Note  from  (2. 7.1-6)  that  V™  =  ^ ns  ■) 

A  suitable  recurrence  relation  is  [Cefola  and  Broucke,  1975] 


Vr 


n+2,s 


(n  —  s  +  1) 
(n  +  s  +  2) 


Vr 


n,s 


(i) 
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Appropriate  initialization  is  provided  by 


bo,o  —  1 
14+1, s+1 


2s  +  2 


K 


(2) 


2.8.3  Calculation  of  Qns  Polynomials 

The  Qns  polynomials  defined  by  (2.8. 1-4)  are  derivatives  of  the  Legendre  polynomials  eval¬ 
uated  at  7,  which  from  (2.2-7c)  is  the  cosine  of  the  angle  between  R3  and  the  angular 
momentum  vector  of  the  satellite.  The  polynomial  Qns  can  also  be  expressed  in  terms  of  the 
associated  Legendre  function  Pns : 


Qns{l)  =  (1  ~72)  s/2Pns( 7) 


(1) 


(Note  from  (2. 7.1-6)  that  V*s  =  Qns( 0)  and  from  (2. 7.4-1)  that  P°s_s{ 7)  =  2s -n^yQnsi'l)  •) 
Recursion  formulas  for  the  Qns  polynomials  follow  directly  from  standard  recursion  for¬ 
mulas  for  the  Pns  functions  [Cefola  and  Broucke,  1975]: 


Qn,s{  7)  1 


f  (2s  -*  l)<5s— i  (7) 

(2s  +  1)7QSjS(7) 

(2 n  -  l)7Qn-i,s(7)  -  (n  +  s  -  l)Qn_2,s(7) 


(n  —  s) 

These  recursion  formulas  are  simply  initialized  by 

Qo,o  —  1 


for  n  —  s 
for  n  —  s  +  1 

for  n  >  s  +  1 


(2) 


(3) 


3  First-Order  Mean  Element  Rates 

As  we  have  seen,  the  first-order  mean  element  rates  Aia  are  given  by  equations  (2.4-18).  The 
osculating  rate  functions  Fia  for  a  conservative  perturbation  are  given  by  (2.2-10)  and  for  a 
nonconservative  perturbation  by  (2.2-5).  In  this  chapter  we  record  the  specific  form  of  these 
equations  for  each  of  several  perturbations. 

3.1  Central-Body  Gravitational  Zonal  Harmonics 

For  the  central-body  gravitational  zonal  harmonics,  the  appropriate  averaging  operator 
<  •  •  •  >  is  (2.4-10)  and  the  appropriate  disturbing  function  7 Z  is  (2.7.1-16)  with  m  =  0. 
Further  details  of  the  reduction  of  the  averaged  equations  to  the  form  recorded  here  may  be 
found  in  [Cefola  and  Broucke,  1975]  and  [McClain,  1978]. 

The  first-order  contribution  of  the  central-body  gravitational  zonal  harmonics  to  the 
averaged  equations  of  motion  (2.4-1)  is 
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da 

dt 

dh 

dt 

dk 

dt 

dp 

dt 

dq 

dt 

dX 

dt 


0 


BdU  k 
~A~dk  +  AB 
BdU  h 
~A~dh 
C 


AB 


ipU,cn  IqUjP 7 ) 

(pG^a 7  -y 


2  AB 
IC 


Uif} 7 

2ABU'ai 
2a  dU 
A  da 


B  dU  ,dUx 
A(i  +  B){h~dh  +  kdk]  + 


AB 


(jPU IOC 7  lQ.Uif3  ry  ) 


(1) 


Here  (a,h,k,p,q,  A)  are  now  the  mean  elements  and  U  is  the  mean  disturbing  function 
(2. 5. 5-2).  In  deriving  (1)  from  (2.2-10),  we  have  made  use  of  the  following  property  of  the 
cross- derivatives  for  the  mean  disturbing  function: 


Ujhk  U,ap  0 


(2) 


The  mean  disturbing  function  reduces  to 


U 


N—2  N 


£  £  (2  -  So.) 

s= 0  n=s+2 


d n.  I  7 


ns 


— n—  l,s 


Qns^d  s 


Here 


Jn  —  Cn  0 

Vns 


K , 


—n—l,s 


Qns(x) 

r<  _ 

^ s  ~  Ltq s 


Since 


geopotential  coefficients 

coefficients  calculated  from  (2. 8. 2-1, 2) 

kernels  of  Hansen  coefficients  calculated  from  (2. 7. 3-6) 

functions  calculated  from  (2. 8. 3-2, 3) 

polynomials  calculated  in  Section  2.7.5 


(3) 


Gs  +iHs  —  Gqs  +  ?77qs  —  {k  +  ih)s  (o:  —  iP)s 


[C9(k,  h)  +iSs(k,  h)][Ca{a,P)-iS9{a,P)]  (4) 


an  alternate  set  of  recursion  formulas  for  the  Gs  polynomials  is 

Gs  =  (ha  +  h/3)Gs-i  —  (ha  —  kP)Hs_  1  ,  Go  =  1 
Hs  =  (ha  -  k(3)Gs_  1  +  (ka  +  hp)Hs_l  ,  H0  =  0 

Note  that  the  Gs  polynomials  are  of  degree  s  in  the  eccentricity;  for  small  eccentricity  orbits, 
the  series  (3)  may  be  truncated  by  prescribing  the  maximum  possible  value  of  the  D’Alembert 
characteristic  s. 
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In  equations  (1),  we  need  the  partial  derivatives  of  U  with  respect  to  (a,  h,  k,  a,  (3, 7). 
These  are  easily  obtained  by  differentiating  equation  (3): 

^7  =  |iE(2-<W(n  +  l)(£f  JnVnsK^sQnsGs 


—  =  —  —  E(2  —  <50s)  (  —  )  JnVnsQns  (  K( 


— n— 1,5 


V7  =  -^E(2-50s)  -)  JnVnsQns  l  Kt 
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«  V  a  ) 

E(2  -  (-) n  JnVnsIun~hsQ 

a  yr  V  a  / 


£(2  -  <So 


JnMns  -^0 


-n-l,sn 

0  9a 

-n-l,sn  d)Gs 
0  Wn.s  ^ 

— n— l,s  dQns 
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Here  we  have  obtained  the  partial  derivatives  with  respect  to  h  and  k  of  K0  n  1,s (y)  from 
the  chain  rule  and  the  definition  (2. 7. 3-4). 

Recursion  formulas  for  0  d (x  are  obtained  by  differentiating  the  recursion  formulas 
(2. 7.3-5): 


(l  +  2s)y2* 

2s 

(n  ~  1)X2 

(n  +  s  —  1  )(n  —  s  +  1) 
+  2  M 

X 


d,Kkn’s 


for  n  —  s 
for  n  —  s  +  1 


dATn+1’s 


for  n  >  s  +  1 


Recursion  formulas  for  7  are  obtained  by  differentiating  (2.8. 1-4): 

^(7)  =  Qn, 8+1(7)  (8) 

Recursion  formulas  for  the  partial  derivatives  of  Gs  may  be  obtained  by  differentiating  (4): 

dG°  nr  tt 

oh 
&  g 

— — —  =  s«Gs_i  +  s(3Hs_  1 

GC  (9) 

wr -  =  s&G's-i  -  shHs_i 


—  shGs_\  -f-  skHs_i 
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If  we  retain  only  the  J2  term  in  the  expansion  (3),  the  central-body  gravitational  dis¬ 
turbing  function  simplifies  to 

^(72  - 


U  = 


where 


a3(l  -  h2  -  A;2)3/2 

j  _  3/ii?2  J-2 


The  contribution  of  (10)  to  the  averaged  equation  of  motion  (1)  is 

da 


dt 

dh 

dt 

dk 

dt 

dp 

dt 

dq 

dt 

dX 

dt 


=  0 


Jk[  3y2  —  1  +  2y  (pa  —  Iq/3)\ 
AB4a3 

Jh[  3y2  —  1  +  2y  (pa  —  Iq/3 )] 
AB4a 3 

CJ/?7 
~  AB4  a3 
ICJa'y 
~  AB4a3 

J[(  1  +  B){ 3y2  -  1)  +  2y (pa 


m) 


(10) 


(11) 


AB4a 3 


In  order  to  update  the  orbital  elements  in  a  differential  corrections  procedure,  it  is  nec¬ 
essary  to  compute  the  partial  derivatives  with  respect  to  the  mean  orbital  elements  of  the 
mean  element  rates  (the  ^  in  the  matrix  F  defined  by  (2.6-7)).  The  partial  derivatives  of 
the  J-2  contribution  (11)  to  the  mean  element  rates  are  easily  obtained,  using  (2. 1.6-1)  and 
(2. 1.9-4).  The  nonzero  derivatives  are 
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dh 

da 

dh 

dh 

dh 

dk 

dh 

dp 

dh 

dq 

dk 

da 

dk 

dh 

dk 

dk 

dk 

dp 

dk 

dq 


7Jk[3 72  —  1  +  2'y(pa  —  Iq/3)\ 

2  aAABA 

AJhk[3 72  —  1  +  27  (pet  —  Iq/3)\ 
a3AB6 

J(  1  -  h2  +  3A;2)[372  -  1  +  2 7(pct  -  Jg/3)] 
a3AB6 

2  J/c[6«7  +  2p(ct2  —  72)  —  2t/2ct7  —  2  Iqf3{a  +  £>7)  +  Cct7] 

a3ABAC 

2IJk[Q(3^  +  2pct7  +  21  q^2  —  2Ipqary  —  2p2/3 7  +  C/37] 

a3ABAC 

7J/i[372  —  1  +  27  (pa  —  ig/3)] 

2aAABA 

(1  -  A;2  +  3/z2)  J[372  -  1  +  27(pct  -  Jg/3)] 
a3AB6 

4J/?,/c[372  —  1  +  27  (pet  —  ig/3)] 
a3AB6 

2JAt[6ct7  +  2p(ct2  —  72)  —  2g2ct7  —  2  Iq(3{a  +  P7)  +  Cct7] 

a3ABAC 

2/J/i[6/?7  +  2pct7  +  2/g72  —  2Ipqa 7  —  2p2(3ry  +  C/?7] 

a3ABAC 
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dp 

da 

dp 

dh 

dp 

dk 

dp 

dp 

dp 

dq 

dq 

da 

dq 

dh 

dq 

dk 

dq 

dp 

dq 

dq 

dX 

da 

dX 

dh 

dX 

dk 

dX 

dp 

dX 

dq 


7CJ £7 
2  aAAB4 
ACJhfii 
a3AB6 
ACJk ^7 
a3AB6 

2J[pf3^  +  a(/3  +  Iq'y)] 
a3AB 4 

2IJ[—ICq/3ry  +  f32  —  72  +  pa  7] 
a3ABAC 

7C'/J«7 
2aAABA 
ACIJha'y 
a3AB6 
ACIJka 7 
a3AB6 

2IJ[pary  +  a2  —  72  —  ig/?7] 
a3AB4 

2J[gC/o;7  +  /3(a  —  jry)] 
a3ABAC 

7J[(1  +  5)(372  -  1)  +  27(pa  -  Iq/3)\ 

2aAABA 

Jh[{ 372  -  1)(4  +  55)  +  87 (pa  -  Jg/3)] 
a3AB6 

JA;[(372  -  1)(4  +  SB)  +  8 7(pa  -  /g/3)] 
a3AB6 

2  J[6(l  +  5)a7  +  2p{a2  —  72)  —  2g2«7  —  2Ipqf3 7  +  Ca7] 

a3ABAC 

2/J[6(l  +  5)/?7  +  2p/3(a  -  ]ry)  +  2/07(7  -  pa)  +  C/37] 

a3ABAC 


(12) 
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3.2  Third-Body  Gravitational  Potential 

For  a  third-body  point  mass,  the  appropriate  averaging  operator  is  (2.4-10)  and  the  appro¬ 
priate  disturbing  function  1Z  is  (2.8. 1-6).  Further  details  of  the  reduction  of  the  averaged 
equations  to  the  form  recorded  here  may  be  found  in  [Cefola  and  Broucke,  1975]  and  [Mc¬ 
Clain,  1978], 

The  first-order  contribution  of  the  third-body  gravitational  disturbing  function  to  the 
averaged  equations  of  motion  is  identical  in  form  to  equations  (3.1-1)  for  the  central-body 
zonal  harmonics.  Of  course,  the  direction  cosines  ( a ,  (3, 7)  have  different  interpretations  for 
the  two  perturbations. 

The  mean  disturbing  function  is  now 


N  N 

c  =  f£  E  (2-W 

3  s=0  n=max(2,s) 


VnsK™QnsGs 


(1) 


Here 


“0 

Qns{ 7) 

Gs 


coefficients  calculated  from  (2. 8. 2-1,  2) 

kernels  of  the  Hansen  coefficients  calculated  from  (2. 7.3-7,  8) 
polynomials  calculated  from  (2. 8. 3-2,  3) 
polynomials  which  can  be  calculated  from  (3.1-5) 


The  partial  derivatives  of  U  needed  in  equations  (3.1-1)  are  easily  obtained  by  differentiating 
equation  (1): 

dU_ 
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d  AT5' 
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Recursion  formulas  for 
(2. 7.3-6,  7): 


dKZ 


dx 


are  obtained  by  differentiating  the  recursion  formulas 


dK™ 

dx 


n— 1,5 


2  n  +  1  dK0 
n  +  1  dx 


0  for  «  =  s-l  or  n  =  s 

(n  +  s)(n  —  s )  dK^~2,s 


n(n  +  l)x2  dx 


+ 


2(n  +  s)(n  -  s )  n_2)S 

7  77 — 7  An 
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n(n  +  l)y3 


for  n  >  s 

(3) 


Recursion  formulas  for  the  derivatives  of  Qns  and  Gs  are  given  by  (3.1-8,  9). 

If  we  wish  to  retain  only  the  dominant  terms  in  the  expansion  (1),  we  should  include  at 
least  the  first  two  terms,  since  the  n  —  2  term  vanishes  in  the  ^  and  ^  equations  for  zero 
eccentricity  orbits,  leaving  the  n  —  3  term  dominant  [Collins  and  Cefola,  1978]. 


3.3  Central-Body  Gravitational  Resonant  Tesserals 

For  the  central-body  gravitational  tessera!  harmonics,  the  appropriate  averaging  operator  is 
(2.4-11)  and  the  appropriate  disturbing  function  is  (2.7.1-16)  with  m  ^  0.  Further  details  of 
the  reduction  of  the  averaged  equations  to  the  form  here  may  be  found  in  [Proulx,  McClain, 
Early,  and  Cefola,  1981]  and  [Proulx,  1982], 

The  first-order  contribution  of  the  central-body  gravitational  tessera!  harmonics  to  the 
averaged  equations  of  motion  (2.4-1)  is  identical  to  (2.2-10),  except  that  77  is  replaced  by  U\ 
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Here  (a,h,k,p,q,  A)  are  now  the  mean  elements  and  U  is  the  mean  disturbing  function 
defined  by 


U  =  <  77  >  = 

+Re 
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47T2  J  —  7T  ./- 

1 


77 (a,  h,  k,  p,  q,  A,  6 ,  t)d\  dO 
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C 3,m)eB 


>  —IT  J  —IT 


(2) 

(Remember  that  B  here  is  the  set  of  all  ordered  pairs  (j,  m)  with  the  properties  (2.4-12,  13)). 

The  mean  disturbing  function  is  identical  to  (2.7.1-16),  except  that  now  only  the  resonant 
tesserals  are  included  in  the  summations: 
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In  equations  (1),  we  need  the  partial  derivatives  of  U  with  respect  to  (a,  h,  k,  A,  a,  /3, 7). 
These  are  easily  obtained  by  differentiating  equation  (3): 

=  Re{-S  E  («  + !)  (7)”  +  iHL) 

°a  a  j,m,s,n  V  a  ' 

( Cnm  -  iSnm)exp[i(j X  -  m9)]j 
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OK  Q  j,m,s,n  V  a  ' 
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Here  we  have  obtained  the  partial  derivatives  with  respect  to  h  and  k  of  KqH~1,s (e2)  from 
the  chain  rule  and  the  relation  h2  +  k2  =  e2,  where  e  is  the  orbital  eccentricity. 

dK™s 

Recursion  formulas  for  ,  \  are  obtained  by  differentiating  the  expansion  (2.7.3-10): 
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Recursion  formulas  for 
(2. 7.4-2,  3): 


are  obtained  by  differentiating  the  recursion  formulas 
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Recursion  formulas  for  the  partial  derivatives  of  GJms  and  PtJrns  are  obtained  by  differentiating 
(2. 5. 3-5)  and  (2.7.5-1,  2): 


s  -  h)Cm_Is{a,(3)  -  I(s  -  j)S\s^l{k,h)Sm_Is{a,f3)  for  |s|  <  m 

s  -  j\Cls_jhl(k,  h)C\a-Im\(a,  j3)  -(s-  j)sgn(m  -  s)Sjs_jhl(k,  h)S\s-Im\(a,  /3) 

for  |s|  >  m 

m 

etc.  Formulas  for  c-^-  are  obtained  by  differentiating  (2.7.1-13). 

In  order  to  update  the  geopotential  coefficients  Cnm  and  Snm  in  a  differential  corrections 
procedure,  it  is  necessary  to  compute  the  partial  derivatives  with  respect  to  the  coefficients 
of  the  mean  element  rates  (the  in  the  matrix  F  defined  by  (2.6-7)).  These  are  easily 
obtained  by  partial  differentiating  (1)  and  (4)  with  respect  to  Cnrn  and  Snrn .  Introducing 
the  parameter 

1  if  max(2,  m,  |s|)  <  n  ,  . 

0  otherwise  ' 

we  can  write  the  results  in  the  compact  form 


r)ii  r)h  2ni  00  N  //?\n 

+  =  E  E  Ci(-)  ImVZK, K-n~uPr(GL  +  iHDexelm - mfl)] 

j—_OQ  \  &  / 

C j,m)eB 

(9) 

etc. 

If  we  assume  that  the  orbital  eccentricity  is  zero,  the  averaged  equations  of  motion  for 
the  resonant  tesserals  significantly  simplify.  See  [Collins  and  Cefola,  1978]. 


3.4  Atmospheric  Drag 

For  atmospheric  drag,  the  appropriate  averaging  operator  is  (2.4-10),  and  the  first-order 
mean  element  rates  are  obtained  by  substituting  (2.2-5)  into  (2.4-18).  To  avoid  having  to 
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solve  Kepler’s  equation,  and  to  smooth  the  perturbation  around  perigee,  we  can  convert  the 
integrals  over  the  mean  longitude  A  into  integrals  over  the  eccentric  longitude  F  or  true 
longitude  L  by  use  of  (2. 5. 2-6)  or  (2.5.3-14).  The  first-order  contribution  of  drag  to  the 
averaged  equations  of  motion  can  then  be  written  in  either  of  the  forms 


da,i 

dt 


1  rp2  fr 
2n  J Fi  \a 


'dai  \ 

wVdF 


(la) 


or 


dai  _ _ 1 _  fL2  fr^2  f  da^  ^  ^ 

dt  27Ta/1  -  h 2  -  k 2  Ai  \a)  \  dr  7 


(lb) 


The  quantities  ^  and  are  given  in  terms  of  the  equinoctial  elements  by  equations 
(2. 1.4-1,  6,  7,  8,  9),a(2.1.6-l)r,  and  (2.1. 7-3). 

The  limits  (F\,  F2)  in  (la)  or  (Li,  L2)  in  (lb)  indicate  the  values  of  F  or  L  at  atmosphere 
entry  and  exit.  If  the  satellite  enters  and  leaves  the  atmosphere  at  a  critical  distance  r  from 
the  center  of  the  central  body,  then 


where 


or 


where 


Fi  =  —E  +  u  +  /n 
F-2  =  E  +  to  +  IQ 


E  =  arccos 


L1  =  -f  +  u  +  m 
L2  =  f  T  oj  T  /O 


/  =  arccos 


a(l-e2)  _ 

r _ 

e 


(2a) 


(2b) 


Of  course,  if  the  satellite  remains  totally  within  the  atmosphere,  the  limits  of  integration 
in  (1)  can  be  taken  to  be  (— n,  7 r). 

The  perturbing  acceleration  due  to  atmospheric  drag  is  commonly  modeled  by  the  formula 
[Escobal,  1965]: 


q  =  %rpiv  ~  _  ^  (3) 

Here 

Cd  =  drag  coefficient  of  satellite 

(Assuming  total  specular  reflection, 

Cd  —  2  for  a  sphere, 

Cd  =  4  for  a  flat  plate  perpendicular  to  (v  —  r).) 

A  =  cross  sectional  area  of  satellite 
m  =  mass  of  satellite 
p  =  density  of  atmosphere 
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r  =  ^  =  velocity  of  satellite 
v  =  velocity  of  atmosphere 

If  we  assume  that  the  atmosphere  rotates  with  an  angular  rate  equal  to  the  angular  velocity 
to  of  the  central  body,  then  v  =  oj  x  r.  The  vector  q  is  resolved  along  the  ( x,y,z )  axes  of 
Figures  1  and  2  for  use  in  the  quadratures  (1). 

The  developers  of  SST  have  used  various  density  models  for  the  upper  atmosphere  of  the 
Earth.  One  of  these  is  the  modified  Harris-Priester  atmosphere  (described  in 
[Long,  Capellari,  Velez,  and  Fuchs,  1989]): 


where 
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Here  pmin(H)  and  pmax(R)  denote  the  minimum  and  maximum  densities  at  a  height  H  above 
a  reference  ellipsoid  <yHl  <  H  <  H2),  Hu  H2,  pmilll,  pmaxi,  pmin2,  Pmax2)  N  are  constants 
whose  values  are  available  from  Tables,  and  by  denotes  the  angle  between  the  diurnal  bulge 
and  the  satellite.  If  b  denotes  a  unit  vector  pointing  from  the  Earth’s  center  to  the  diurnal 
bulge,  then 

b  •  r 

cos  (pb  = -  (5) 

r 

The  diurnal  bulge  follows  the  path  of  the  Sun  but,  because  of  the  Earth’s  rotation,  lags  the 
sub-solar  point  by  an  angle  9 &  (approximately  30°  at  2  P.M.  local  time).  We  can  obtain  the 
vector  b  by  rotating  the  vector  z b  (pointing  from  the  Earth’s  center  to  the  sun)  through  an 
angle  9b  about  the  Earth’s  axis  of  rotation.  Letting  R  denote  the  3x3  matrix  whose  elements 
are  direction  cosines  between  the  (x,  y,  z )  axes  and  an  Earth-fixed  set  of  Cartesian  axes,  we 
can  write  the  transformation  law  between  the  (x,  y,  z )  components  of  b  and  z b  as 


bx 

cos  6  b 
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ZBx 
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Here  R7  denotes  the  transpose  of  the  matrix  R  (see  [Danielson,  1991]). 
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3.5  Solar  Radiation  Pressure 

The  general  equations  expressing  the  first-order  contribution  of  solar  radiation  pressure  to  the 
averaged  equations  of  motion  are  formally  identical  to  the  equations  (3.4-1)  for  atmospheric 
drag. 

The  limits  (T\,  F2)  in  (la)  or  (L1;  L2)  in  (lb)  now  indicate  the  values  of  F  or  L  at  shadow 
exit  and  entry.  If  we  assume  the  shadow  is  a  circular  cylinder  in  shape,  the  shadow  exit  and 
entry  longitudes  are  determined  by  the  shadow  equation  (as  explained  in  [Escobal,  1965]  and 
[Cefola,  Long,  and  Holloway,  1974]): 

5  =  0  (1) 

Here 

S  =  l  —  M.{1  +  k  cos  L  +  h  sinL)2  —  (a  cos  L  +  [3  sin  L)2 


M 

a 

P 


Rl 


a2(l  —  h2  —  k 2) 

R.3 

r3 

R,3 


•  f 


R? 


■  g 


where  is  the  central-body  radius,  and  R3  is  the  vector  from  the  center  of  mass  of  the 
central  body  to  the  sun.  To  obtain  the  solutions  to  equation  (1),  the  following  quartic 
equation  must  be  solved: 

.Ao  cos4  L  +  A\  cos3  L  +  A2  cos2  L  +  cos  L  +  A4  =  0  (2) 


where 


A0  =  4  B2+C2 
A 1  =  8B  Mh  +  4C  Mk 
A2  =  -4B2 +  4M2h2  ~2V  C  +  4M2k2 
=  -8B  Mh  -  4VMk 
Aa  =  -4  M2h2  +  V2 
B  =  a/3  +  M.hk 
C  =  a2  -  (32  +  M(k2  -  h2) 

V  =  1  —  /32  —  M(1  +  h2) 


The  real  roots  of  (2)  must  be  sorted  to  eliminate  extraneous  roots  and  to  determine  the  exit 
and  entry  values  of  true  longitude  L.  The  correct  values  of  L  must  satisfy  (1)  as  well  as  the 
condition 


R3  y 

■  -  =  cos  (j)  =  a  cos  L  +  (3  sin  L  <  0 
R3  r 


(3) 
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At  exit  from  shadow 


(4) 


while  at  entry  into  shadow 


dS 

M>0 


dS 

dL<0 


(5) 


Of  course,  if  the  satellite  remains  totally  within  sunlight,  the  limits  of  integration  in  (3.4-1) 
can  be  taken  to  be  (— 7r,  it). 

The  perturbing  acceleration  due  to  solar  radiation  pressure  is  [Cefola,  1982]: 


.  (r  -  R3) 

2m  c  0  |r  —  R3|3 


(6) 


Here 

Cr  =  radiation  pressure  coefficient  of  satellite 

(Assuming  total  specular  reflection, 

Cr  =  2  for  a  spherical  mirror  or  black  body, 

Cr  =  4  for  a  flat  mirror  perpendicular  to  (r  —  R3).) 
A  =  cross  sectional  area  of  satellite 
m  =  mass  of  satellite 
C  =  mean  solar  flux  at  surface  of  sun 
c  =  speed  of  light 
Rq  =  radius  of  sun 

r  —  R3  =  position  vector  from  sun  to  satellite 
If  we  suppose  that  the  satellite  is  always  in  sunlight,  and  that  the  parameter 


T 


CrAC 
2  m  c 


is  constant,  we  can  derive  (6)  from  the  disturbing  function 


n  = 


Use  of  the  expansion  (2.8-2)  then  leads  to 


n 


T 

R~3 


Pn{  cos  4>) 


(7) 

(8) 

(9) 


The  radiation  disturbing  function  (9)  is  identical  to  the  third-body  disturbing  function 
(2.8-3),  except  that  the  factor  /i3  is  replaced  by  — T  and  the  summation  starts  from  n  =  1. 
Hence  we  can  immediately  write  down  the  mean  radiation  disturbing  function  by  analogy 
with  (3.2-1): 


U 


s~r  N 

—y 


(10) 
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If  we  retain  only  the  first  nonzero  term  in  the  expansion  (10),  the  mean  radiation  dis¬ 
turbing  function  simplifies  to 

U  =  V{ka  +  h/3)  (11) 


where 


V  = 


3  T  a 

2 M 


The  contribution  of  (11)  to  the  averaged  equations  of  motion  (3.1-1)  is 

da 


dt 


=  0 


dh  BVa  V&7 
dt 
dk 
dt 


A  AB 
BV(3  Vhy 


A 

dp  CVh'-f 
dt  2AB 
dq  ICVk-i 
dt 
dX 
dt 


+ 


AB 


( kp  —  Ihq ) 

(. kp  —  Ihq ) 


(12) 


2  AB 

(2  +  B)V{ka  +  h/3)  Vy 
A(l  +  B)  AB 


( kp  —  Ihq ) 


4  First-Order  Short-Periodic  Variations 

Knowing  Fourier  series  expansions  for  the  osculating  rate  functions  Fia  or  the  osculating 
disturbing  function  7 Z,  we  can  construct  the  first-order  short-periodic  variations  r)ia  from  the 
results  in  Section  2.5.  In  this  chapter  we  record  the  specific  forms  of  the  expansions  for  each 
of  several  perturbations. 


4.1  Central-Body  Gravitational  Zonal  Harmonics 

For  the  central-body  gravitational  zonal  harmonics,  the  appropriate  disturbing  function  77 
is  (2.7.1-16)  with  m  =  0.  From  the  results  in  Sections  2.5.1  or  2.5.5,  we  can  construct  a 
Fourier  series  expansion  in  the  mean  longitude  A  for  the  first-order  short-periodic  variations 
Tjia,  as  was  done  by  Green  [1979]. 

However,  for  the  central-body  zonal  harmonics,  it  is  possible  to  construct  a  finite  modified 
Fourier  series  in  the  true  longitude  L  for  the  r)lo .  For  single-averaged  perturbing  forces  which 
increase  rapidly  as  the  satellite  approaches  the  central  body,  notably  central-body  zonal 
harmonics  and  atmospheric  drag,  a  Fourier  series  expansion  in  L  will  converge  faster  than 
equivalent  expansions  in  F  or  A  when  the  eccentricity  of  the  satellite  orbit  is  large.  This 
happens  because  the  magnitude  of  such  a  perturbation  is  strongly  peaked  around  perigee 
when  considered  as  a  function  of  A.  Transforming  the  independent  variable  to  L  broadens 
the  peak  considerably,  making  it  easier  to  approximate  with  a  finite  Fourier  series.  In  this 
Section  we  outline  the  construction  of  this  most  desirable  expansion  in  L.  Further  details  may 
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be  found  in  [Cefola  and  McClain,  1978],  [Kaniecki,  1979],  [McClain  and  Slutsky,  1980],  and 
[Slutsky,  1980]. 

The  disturbing  function  less  it’s  mean  can  be  written  as 


Here 


{..  N  n  ,  td \  n 

--EE(2-5  os)  -  JnVnsQns[Cs(a:P)  +  iSs(a,P)} 

a  n= 2 s=0  '  °  ' 


E  Yj  ”  ’  "exp (ijX) 


J=-oo 
3+ 0 


n  0 


Jn  —  —Ci 

Vns 

Qns('y) 

Cs(a,/3),  Ss(a,P) 

yr—n—lj  —  S 


geopotential  coefficients 
coefficients  calculated  from  (2. 8. 2-1, 2) 
polynomials  calculated  from  (2. 8. 3-2,  3) 
polynomials  calculated  from  (2. 5. 3-6) 
coefficients  of  the  expansion  (2. 7. 1-8): 


(1) 


/  n  \  n+ 1  °o 

f-J  exp (-isL)  =  yj~n~h~sexp(ijA) 

' r '  j=- oo 


(2) 


The  short-periodic  generating  function  (2. 5. 5-4,  5)  is  easily  obtained  by  integrating  (1): 


!.  N  n  /  R\n  00  Y~n  ^’~Sf 

--EE(2-^0s)  -  JnVnaQns[Cs(a,P)  +  iSs(a,P)]  E  ~ ^ 

a  n= 2 s=0  '  a '  j=-o o  C 

3+ 0 

(3) 

The  infinite  series  in  the  mean  longitude  A  in  (3)  can  be  replaced  by  a  finite  modified 
series  in  the  true  longitude  L.  To  see  this,  first  integrate  both  sides  of  (2)  with  respect  to  A 
to  obtain 


oo  Yrn~1’~sexp(ij\) 


E  ^ 

j=-oo 
3+ 0 


*7 


rA  /  n\  n+l 

-  )  exp(— isL)d\  —  AY'C 
r  ) 


— n—  1,— s 


Next  perform  the  integral  in  (4)  by  using  the  expansion 


=  \/l  —  h2  —  k2  ^2  Yo  “  ’  Jexp  (ijL) 


3=—n 


and  the  change  of  variable  (from  2.5.3-14) 


dX  = 


VI  -  h2  -  k2  \aj 


r\ 


-  dL 


(4) 

(5) 

(6) 
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The  infinite  series  then  becomes 


~  Yj  n  h  sexp(ij\) 


"-1  yo~n-1,-Jexp[i(j 

=^-i)  i(<j  ~  s) 


s)L] 


+  y0"n"1,_s(L  -  a) 


(7) 


Replacing  the  infinite  series  in  (3)  with  (7),  and  introducing  the  kernels  K0  n  lj  of  the 
Hansen  coefficients  through  (2.7.1-11),  we  obtain 

{.  N  n  /  R  \  n 

-EE(2-W  -  JnVnsQns 

O'  n=2  s= 0  \  0/  / 

Y  [C8{a,P)  +  iS8(a,P)][C\j\(k,  h)  -  isgn(j)S\j\(k,  /r)]/l0~ra~ljexp[i(j  -  s)L\  1  (8) 

j-Tt-i)  “  s>  J 

i+s 


Since  all  of  the  symbols  in  (8)  except  for  i  =  \/— T  are  real,  we  can  easily  cast  S'  into  the  real 
form 


s  =  u(l  —  a)  —  ^  y  y  (2  -  *»)  ;  y  x, 


N 


n—  1 


-n-lj 

0 


71—2  s=0 


,l=o 

It^s 


cos(j  —  s)L  +  GjS  sin(j  —  s)L 
j  ~  s 


where 


Kin) 

Gis 


H 


JS 


J 


11—1 

E  Kon~1J 

1=1 


-Ijs  COS  (j  +  s)L  +  JjS  sin  (j  +  s)L 


j  +  s 


(9) 


/  D\  n 

J  VnsQnsil) 

Cj(k ,  h)Cs(a ,  /?)  +  S^/c,  h)Ss(a ,  /3) 

Cj(k,h)Sa(a,/3)  -  Sj(k,h)Cs(a,l3)  (10) 

Cj(k ,  h)Ss(a ,  /3)  +  S'j(/c,  h)Cs(a ,  /?) 

Cj(k,h)Cs(a,l3)  -  Sj(k,h)Ss(a,  (5) 


The  last  step  is  a  redefinition  of  indices,  so  as  to  isolate  the  coefficients  of  cos  jL  and 
sin  jL,  and  a  rearrangement  of  the  order  of  summation.  The  short-periodic  generating  func¬ 
tion  can  then  finally  be  written  as  the  finite  modified  Fourier  series 


S 


2N-1 


C°  +  U{L-  A)  +  E  K  cos  JL  +  &  sin  JL ) 

l=i 


(11) 


Here  the  coefficient  C°  is 


C° 


2N-1 

y  (CVj+SVy) 


1=1 


where  pj  and  a3  are  given  by  (2. 5. 3-3,  4).  The  coefficients  &  are 


(12) 
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N- 1  N 


C>  =  -fr  E  E  ( 

J  y  s=j  n=s+l 

N-j  TV  ] 

-  E  E  (2  -  S„j+,)JnH,j+tK^-1-’Vn+-  -  2X?U)JjHt,ijK^-'-0Ljl 

s= 0  n=max(j+stj+l) 


+77  {  Z?U)  E(2  -  s« J-.W.J-.K^L 

aJ  S=1 


-3- 1,«  r  js 

o  -k? 


j  N 


+iE‘(i)  E  E  (2  -  +x|'v-1(j) Ei 

s=l  n=j+l 


where 


2  1  min(j-l,W) 

Ei  =  E  E  (2  -  S0J_,)JnI,ij-,K^1-‘Li- 

s=j—min(j—l,N)  n=j—s 
minQ— 1,7V)— 1  min(j  —  1,7V) 

+  E  E  (2-  5o,j-s)JJs,j-sKQ  n^1'sL3~s  for  j  even 

n=s+l 


(14a) 


minQ— 1,7V)— 1  minQ— 1,7V) 

Ei  =  E  E  (2  -  S0J-,) 

- 7  —  1  n=S+l 

*  2 

f  0  for  N  =  2,3 


LT'  min(j  — l,iV) 


for  j  odd 


E  E  (2  -  Soj-sWsj-sKo^Li-8  for  AT  >  4 

s=j— min(j— 1,N)  n=j—s 


Similarly,  the  coefficients  S 3  are 


S3  =  - 


N- 1  AT 


-IxEti)  E  E  (2  -  6o,,-j)JnG,,,-iK^n-1’‘L--i 

0*3  y  s=j  n=s+l 


TV— j  TV  1 

E  E  (2  -Saj+,)JnG,J+,K^-l-‘Li+‘  +2lZ(j)J1G„jK^-l’°Li 

s=0  n=max(j+s  j+1) 
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-4  xf(i) 

aj 


j'-i 


E(2  -  50j_s)  J,- Jsj-sXo"j'-1,s^“s 


+*1  tt) 


S=1 

j  N 


E  E  (2-50j-s)JnJSJ-S^o'n_M"“S 


5  TJ~ 
^ n 


s=l  71=7+ 1 


where 


Ej2 


2_1 

2 


— n—  l,s  T  j—s 

0  L'n 


min(j  — 1,7V) 

E  E  (2-  <S0 ,j-s)JnJs,j-sK( 

s=j—mm(j—l,N)  n=j—s 
min(j— 1,7V)— 1  min  (j— 1,7V) 

E  E  (2  -  Soj-JJnJsj-sK^-^LL 

n=s+ 1 


s=2 
5  2 


for  j  even 


(16a) 


min(j  —  1,7V)— 1  min(j  — 1,7V) 

'3  -  £  £  (2  -  S0j.,)JnJ,j.,K^'-‘Vn- 

j- 1  n=s+l 

6  2 


EJ2  = 


(  0  for  N  =  2,3 


+ 


for  j  odd 


X 


-27V— .3 


^  min  (j  — 1,7V) 


(7)  E  E  (2  -  8Qj-s)JnJ8j-sK{ 

s=j—mm(j—l,N)  n=j—s 


-n-i,,Lj-.5for  JV  >  4 


(166) 

Note  that  the  hrst  index  of  the  G,  H,  /,  J  polynomials  defined  in  (10)  indicates  their  degree 
in  the  eccentricity;  for  small  eccentricity  orbits,  the  series  (14)-(16)  may  be  truncated  by 
prescribing  the  maximum  possible  value  of  s. 

The  first-order  short-periodic  variations  r]ia  generated  by  the  function  S  given  by  (11) 
can  be  derived  using  equations  (2.2-10),  (2. 5. 2-7),  (2.5.5-10),  and  the  following: 


dL 

~dh, 


dL 

dk 


Er3 

B3 


kb  +  kb  ^1  +  ^  ^  +  2  (B  +  k2b )  cos  L 


k  h  'i 

+26X6  sin  X  +  —  [B  +  {k2  —  h2)b\  cos  2 L  +  —  (X?  +  2k2b)  sin2X) 


Ef3 

B3 


{h,B  +  hb  ^1  + 


h2  +  k2' 


+  2hkb cos  L 


(17) 


+2 (X?  +  h2b)  sinX  +  ^[— B  +  ( k 2  —  h2)b\  cos  2 X  +  ^-(X?  +  2626)  sin2X) 

r)T  1  ,9  x  /i2  x  P  h2  —  h2 

—  =  — (  +2fccosX  +  2fesinX  +  - - -  cos2X  +  6Xsin2X| 

d\  B*  t  9  9  J 


Here,  as  usual,  h  and  k  are  equinoctial  orbital  elements  and  6  and  B  are  defined  by 
(2. 1.4-4)  and  (2.1.6-lb). 
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In  the  absence  of  explicit  time-dependence,  the  first-order  short-periodic  variations  can 
be  written  as  the  finite  modified  Fourier  series 

2N+1 

Via  =  C°  +  Di(L  -  A)  +  53  ( C{  cos  jL  +  Si  sin  jL)  (18) 

i=i 

(Remember  that  the  equation  of  the  center  (L  —  A)  may  be  calculated  from  (2. 5. 3-2, 3, 4).) 
Expressions  for  the  coefficients  in  (18)  are  given  below  in  terms  of  the  following  quantities: 

1.  The  equinoctial  elements  (a,  h,  k,p,  q )  and  the  retrograde  factor  /  (equation  (2. 1.2-2)). 

2.  The  direction  cosines  (cc,  /?, 7)  of  the  perturbation  symmetry  axis  in  the  equinoctial 
reference  frame  (equations  (2. 1.9-1)). 

3.  The  Kepler  mean  motion  n  (equation  (2. 1.4-3)). 

4.  The  auxiliary  parameter  x  (equation  (2.7.3-4)). 

5.  The  mean  disturbing  function  U  of  the  perturbation  (equation  (3.1-3)). 

6.  The  coefficients  and  of  the  modified  Fourier  series  expansion  in  L  for  the  short- 
periodic  generating  function  S  (equations  (13)-(20)). 

7.  The  cross-derivative  operator  (equation  (2.2-8)). 

8.  The  inclusion  operator  (equation  (2.5.3-18)). 

The  constant  terms  in  the  expansions  (18)  are  given  by: 

2N+1 

c?  =  -  E  (C’iPi  +  Siaj)  (19) 

3= 1 
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The  short-periodic  coefficients  for  the  semimajor  axis  a  are  given  by: 


C{  = 
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The  short-periodic  coefficients  for  element  h  are  given  by: 
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The  short-periodic  coefficients  for  element  k  are  given  by: 
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( kC1  +  hS1 


hU 


-AN~3U) 

-AN~2U) 


X 


'-An2  a2 
X 


(. j  +  2)(kCj+2  +  hSj+ 2) 


n2a 2 

r2N—  1  /  •  \  f  ^X 


-2r_1(j) 


Ln2a2 


(j  +  1)C'+1 


(P5,y 


+ 


3/cy  ■ 

+ 


131 J  '  yn2a2  dh  '  2n2a2 

Ln/Gr  J 

X  0'-2)(-lfeCJ'-2  +  /i5J'-2) 


+^JV+1(i) 


,0.  = 


1  dU 


An2  a2 

hx 


(pU,a 7  lQ.U,p^  ) 


Xn2a2  dh  n2a2Vr~,Q7 
The  short-periodic  coefficients  for  element  p  are  given  by: 


CZ  =  2r_1(j) 


2  n2a2 


-c,y  +?(ci  -c,i„  -jsj 


s{  =  xp-'u) 


2jv-i/ ■ -7  (!  +P2  +  <?2)xr 


2n2a2 


~S,L  +P(sik  -SS 


aft 


da  = 


(1  +  p2  +  q2)xU,(3y 


2  n2a2 

The  short-periodic  coefficients  for  element  q  are  given  by: 


ci  =  AN~\j) 


2n2a2 


-ICL 


i,+q(cik-Ci^3Si) 


Si  =  AN-\j) 


-i/ (f  +P2  +  <?2)x  r 


2n2a2 


- IS,2n+q(S,ik-S,^+jC ■>) 


D,  = 


(l+p2  +  q2)xIU, 


o.'y 


2  n2a2 


(22) 


(23) 


(24) 
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The  short-periodic  coefficients  for  the  mean  longitude  A  are  given  by: 


ci  = 


-AU) 

-nu) 


X 


L 2n2a2(l  +  x) 
2 


(4  hU  +  k '2  1rCl  +  hkS1) 


X 


-hkU 


'-n2a2(  1  +  x) 


-jz 

xi 


2N—2  1 


Wa2(l  +  x)y  2 

x\]  +  1) 


^Mrri)(kCi+1+hSi+1) 


+Tf-\j) 


2  dCj 


+ 


X 


n2a  da  n2a2(  1  +  x)  ^  9h  dk  1  n2a2 


n2a2 


Cj 


Si  = 


+zT+1(i) 

i!(i) 


x2(i-2) 


C>-‘  -  hkS> 


-2 


2n2a2(l  +  x)V  2 

V2  /  t2  —  /l2 

1  -  /i/tc1  +  - —  S1 


L2-n2a2(l  +  x) 


(25) 


+Zf(j) 


r  X 

k  -  h 

-n2a2(l  +  x) 

2 

[/ 


+xr~3(j) 


X2(J  +  2)  (hkCj+2  -  k—JlS1+2) 


-2n2a2(l  +  x) 


+X2N~1(j) 


2  1  / ,  dSj  .  dSj  N  \  \  3 


+ 


n2a  da  n2a2(  1  +  x)  ^  d/i  <9fc  '  n2a 


(hm7  +  kmr)  +  ^sX~^X) 


n2a 2 


SJ 


+I“w[sSfrrt)(,lC,~,+A5i") 

X2(J  -  2) 


D,  = 


+llN+\j) 

2  dU 


(hkCj~2  +  ^Z^-Sj-2 


L2n2a2(l  +  x)v .  2 

1  /.  dU  ,  dU \  X 


n2a  da  n2a2 


(TTx)(h~dK  +  k~dk> +  ^(pU^  ~,qU^) 
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Note  that  the  D,,  coefficients  are  simply  related  to  the  first-order  mean  element  rates  Aia 
(given  by  the  right  sides  of  (3.1-1)): 

Di  =  —  (26) 

n 

From  the  central-body  gravitational  potential,  there  are  three  possible  sources  of  explicit 
time-dependence: 

1.  Motion  of  the  central-body  symmetry  axis.  For  the  Earth,  this  is  a  combination  of 
precession  of  the  equinoxes,  nutation,  and  polar  motion. 

2.  Variations  in  the  central-body  rotation  rate. 

3.  Tidal  potential. 

The  principle  effects  of  1  and  2  are  accounted  for  in  SST  by  using  at  each  time  step  the  epoch 
triad  (xB,yB,zg)  to  evaluate  the  direction  cosines  (a,  (3, 7)  from  (2. 1.9-1)  and  the  rotation 
angle  6  from  (2. 7. 1-3, 4).  However,  for  the  Earth  the  above  sources  are  thought  to  be  too 
small  or  two  slowly  varying  to  cause  significant  explicit  time-dependence  effects. 

In  order  to  update  the  orbital  elements  in  a  differential  corrections  procedure,  it  is  nec¬ 
essary  to  compute  the  partial  derivatives  with  respect  to  the  mean  elements  of  the  short- 
periodic  variations  (the  in  (2.6-3,  4)).  The  partial  derivatives  of  the  Jj  contribution  to 
(18)  are  currently  available  in  the  SST  code,  for  the  special  case  of  zero  eccentricity  and 
replacement  of  (a,  /3, 7)  with  the  explicit  formulas  (2. 1.9-3)  in  p  and  q  (thus  motion  of  the 
central-body  symmetry  axis  is  neglected). 

4.2  Third-Body  Gravitational  Potential 

For  a  third-body  point  mass,  the  appropriate  disturbing  function  1Z  is  (2.8. 1-6).  From  the 
results  in  Sections  (2.5.1)  or  (2.5.5),  we  can  construct  a  Fourier  series  expansion  in  the  mean 
longitude  A  for  the  first-order  short-periodic  variations  r)ia,  as  was  done  by  Green  [1979]. 

However,  for  the  third-body  disturbing  function,  it  is  possible  to  construct  a  finite  Fourier 
series  in  the  eccentric  longitude  F  for  the  r)ia .  Since  the  D’Alembert  characteristics  are 
bounded,  this  solution  is  of  closed  form  in  the  eccentricity.  In  this  Section  we  outline 
the  construction  of  this  most  desirable  expansion  in  F.  Further  details  can  be  found  in 
[McClain,  1978],  [Cefola  and  McClain,  1978],  [Slutsky  and  McClain,  1981],  and  [Slutsky, 
1983], 

The  disturbing  function  less  its  mean  can  be  written  as 

{..  N  n  ,  \n  00 

^EE(2-<W  by  K sQns[Cs(a,(3)  -iSs(a,P)]  £  Y™exp(ij\) 

K3n=2s= 0  j=- oo 

3+  0 

(1) 

Here 
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Vns  =  coefficients  calculated  from  (2. 8. 2-1,  2) 

Qns( 7)  =  polynomials  calculated  from  (2. 8. 3-2,  3) 

Cs(a,P),  Ss(a,/3 )  =  polynomials  calculated  from  (2. 5. 3-6) 

=  coefficients  of  the  expansion  (2. 7. 1-8) 

The  short-periodic  generating  function  (2. 5. 5-4,  5)  is  easily  obtained  by  integrating  (1): 

N  n  '  -  x  n  00  yn-sexp(i)A) 


S  =  Rei^HH(2-'5»*)(-5-)  Vn,Qn,{C,(a,P)-iS,(a,P)]  £ 


R 


3  n= 2 s=0 


R, 


j=-oo 
3+ 0 


*7 


(2) 


The  infinite  series  in  the  mean  longitude  A  in  (2)  can  be  replaced  by  a  finite  series  in  the 
eccentric  longitude  F .  To  see  this,  first  integrate  both  sides  of  (2. 7. 1-8)  to  obtain 


g  F/sexp(ijA)  /-a  (r 


3=— 00 
jVo 


*7 


/A  /  rr*  \  TL 

(  —j  exp(isL)dX  —  XYq 


Next  perform  the  integral  in  (3)  by  using  the  expansion 


(3) 


-j  exp  (isL)  =  ^2  WTsexp  (ijF) 

a  j=-n 


where 


w™  =  y0n_1,s 

and  by  using  the  change  of  variable  (from  2. 5. 2-6)) 


(4) 


(5) 


dX  =  (  -  )  dF 

My 


The  infinite  series  then  becomes 


^  P“exp(ijA)  W;+,’%xp(ijF)  ,  „ 

+yo  A) 

j=  00  j=-(n+ 1) 

jyo  iyo 


(6) 


(7) 


Replacing  the  infinite  series  in  (2)  with  (7),  recalling  the  relationships  (2. 5. 3-5),  (2.7.1-11), 
(3.1-4),  and  introducing  the  mean  disturbing  function  U  through  (3.2-1),  we  obtain 


5  =  U(F-  A) 


|T,  )  /'3  ^  t  n  [^(oj,/?)  —  *5,s(af,/3)]W7+1’®exp(*jJF) 

+Re  \  Z.  ZJ2  “  )  K.sQn.s  Z^  - —  -  - 


3  n= 2 s=0 


R;- 


j=-(n+ 1) 


U 


(8) 
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(9) 


The  coefficients  W™s  may  be  expressed  in  the  form 

w;s  =  wY[C\j-s\{k,  h )  -  isgn(j  -  s)S^sl(k,  h)] 

The  functions  WjS(e)  possess  the  Jacobi  polynomial  representation 

(n  +  s)!(n  —  g)!  2\-|s| p\j-s\,\j+s\,  a  f  ,  ,  >  , 

(//-  ./)!(//  -  j)\  C)  ^  {X)  t0r|s|-|j| 

(1  -  c*yWP^f3+s\x)  for  |s|  <  \j\ 


wr(e)  =  []—AA  (-c)li_s|  <!  (n  +  j)\(n  -  j)\K±  1  n  |s| 


1  +  c2 


Here  again  x  is  defined  by  (2. 7.3-4),  e  is  the  orbital  eccentricity,  and 


c  = 


Vh2  +  k2 


1  +  a/1  —  e2  1  +  a/1  —  h2  —  k2 


(10) 


(11) 


After  substituting  (9)  into  (8),  we  can  easily  cast  S  into  its  real  form.  The  indices  may 
be  rearranged  as  usual,  and  Kepler’s  equation  (2.1.4)  may  be  used  to  replace  (F  —  A)  in  (8) 
by 

F  —  X  —  ksinF  —  hcosF  (12) 

The  short-periodic  generating  function  can  then  finally  be  written  as  the  finite  Fourier  series 


jV+l 

S  —  C°  +  U(k  sin  F  —  h  cos  F)  +  ^  (Cj  cos  jF  +  Sj  sin  jF) 

i=i 

Here  the  coefficient  C°  is  (implied  by  <  S'  >=  0  and  (2. 5. 2-8)) 

o  k  i  h  , 


(13) 


Cu  =  -C1  +  -5 
2  2 


(14) 


The  coefficients  and  SJ  are 

N  N  q 


cj  =  J2  J2  —  {-e  b  s]Wj+1,s[sgn(j  -  s)Cs(a,(3)S]:j_s\(k,  h)  +  Ss{ol, (3)C\j_s\(k,  h)\ 


s=0  n=max(2j— l,s)  7 

+e~^wlY-‘[-C,(a,  0)Sj+,(k,  h)  +  S,(a,  p)Cj+,(k,  ft)]} 


(15) 


N  N 


Gr 


s1  =  —p{e  |j  slw]+1,s[Cs(a,(3)C\j_s\(k,  h)  -  sgn(j  -  s)Ss(a,  (3)S]:j_s](k,  h)] 

(16) 

Gns  =  ^(2  -  Jos)  (^)n  VnsQns  (17) 


s=0  n=max(2j  —  l,s)  7 

+e-^w^u[Cs(a,(3)Cj+s(k,h)  +  Ss(a,(3)SJ+s(k,h)}} 


where 
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The  first-order  short-periodic  variations  rjta  generated  by  the  function  S  given  by  (13)  can 
be  derived  using  equations  (2. 1.8-3),  (2.5.5-10),  and  the  following  (obtained  by  differentiating 


(2. 1.4-2)): 

d  F 

~dh 

dF 

dk 

dF 

~d\ 


a  T? 

—  cos  t 
r 

—  sin  F  (18) 

r 

a 

r 


where  r  is  given  by  (2. 1.4-6).  The  partial  derivatives  of  S  with  respect  to  the  elements 
a,a,[3, 7  follow  by  straightforward  differentiation  of  (13): 


dS 


dC° 


dU 


d(a,a,(3 ,7)  d(a,  cc,  /3, 7)  d(a,  a,  (3, 7) 


( k  sin  F  —  hcosF ) 


/  dCj 
\<9(a,  01,  P,  7) 


cos  j  F  + 


dS3 


d(a,a,P,  7) 


sin  jF 


(19) 


The  coefficients  C3 ,  S3  and  the  mean  disturbing  function  U  are  functions  of  the  semimajor 
axis  a  through  the  powers  of  the  parallax  factor  alone,  functions  of  the  direction  cosines  a 
and  (3  through  the  polynomials  Cs(a,fi),  Ss(a,/3),  and  functions  of  the  direction  cosine 
7  through  the  polynomials  Qns(l)-  A  finite  Fourier  series  representation  for  may  be 
obtained  by  partial  differentiation  of  (2)  with  the  infinite  series  replaced  by  (2. 7. 1-8),  followed 
by  substitution  of  (4)  and  (9)  and  the  usual  reduction  to  the  real  form 

oc  N 

=  C°\  +  E(C?a  cos  jF  +  SfA  sin  jF)  (20) 


Here  the  coefficients  are 


,a 


& ’x 


k  h  , 

2C’a  +  25a 
N  N 


E  E  G'n4e-|^s|Wf[C's(a,/3)qj_^A;,/i)^sgn(j-S)^(a,/3)(Sb-_s|(A;,/i)] 


s— 0  n=max(2  ,j,s) 


six 


+e-ti+*)w™ \[Cs(a,  (3)Cj+s(k,  h)  +  S„(a,  f3)SJ+s(k,h)}} 


(21) 


N  N 


E  E  Gns  {e-V-Uw]3 [sgn(j  -  s)C>,  frS^k,  h)  +  S8(a,  /3)C\j.s\ (k,  h)} 


s— 0  n=max(25i7,s) 


+e-^w™[Cs(a,f3)Sj+s(k,h)  -  Ss(a,  (3)Cj+s(k,  h)]} 


(22) 
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The  partial  derivatives  of  S  with  respect  to  the  elements  h  and  k  may  be  obtained  by 
differentiation  of  (13)  and  use  of  (18)  and  (20): 


dS  dC°  dU  .  ^ 
M  =  dh+kWSmF~ 


N+l 

E 

3= 1 


'da 

dh 


d(hU) 

dh 


+  C°A 


cos  F 


(23) 


1  ,  ,  1  ,  ,  \  fdS> 

a  -a  cosjF+br 


-5i_1 
2  ’A 


-  25?1 )  sin  jF 


dS_ 

dk 


dC° 

dk 


N+l 

E 

3= 1 


d{kU) 

dk 

da 


a 


sin  F 


,  dU 

h-r-  COS  F 

dk 


dk 


~  1S3\  1  +  W1)  cos  jF  + 


'dSi 

dk 


(24) 


3~  1 
,A 


-CfA+1  sinjF 


In  the  summations  in  (23)-(24)  CJA  and  S\  are  defined  to  be  zero  for  values  of  j  outside  the 
range  1  <  j  <  N.  The  dependence  of  C-7,  S J  and  U  on  the  elements  h  and  k  is  through 
the  polynomials  Ce(k,h),  Se(k,h),  the  eccentricity  e,  and  the  coefficients  wY,Kqs.  In  the 
absence  of  explicit  time-dependence,  the  first-order  short-periodic  variations  can  thereby  be 
written  as  the  finite  Fourier  series 

AT+l 

Via  =  c;°  +  £  ( Ci  cos  jF  +  Si  sin  jF)  (25) 

3= 1 

where  the  coefficients  Cf  are  given  by  (2.5.2-15a). 

Green  [1979]  studied  the  effects  of  explicit  time-dependence  by  including  several  terms 
in  the  formulas  (2.5.1-15)  for  the  coefficients  Cj  and  Sj  of  the  A-expansions  (2.5.1-14)  of 
the  short-periodic  variations.  He  used  finite  difference  formulas  to  compute  the  partial 
derivatives  with  respect  to  time  of  the  coefficients  Cj  and  Sj  in  the  A-expansions  (2.5.1-11) 
of  the  osculating  rate  functions.  For  medium  or  high  altitude  Earth  satellites,  he  found 
that  the  Lunar  point  mass  perturbation  varies  quickly  enough  for  explicit  time-dependence 
effects  to  be  significant.  McClain  and  Slutsky  [1988]  also  found  that  the  inclusion  of  explicit 
time-dependence  effects  due  to  the  Moon  and  Sun  improved  the  performance  of  SST  for  high 
altitude  Earth  satellites. 

We  can  include  the  effects  of  explicit  time-dependence  in  the  F-expansions  of  the  present 
Section  by  using  (2.5.1-10)  or  (2.5.5-11)  and  expressions  in  Section  2. 5. 2.1,  the  first-order 
short-periodic  variations  r]f'a  including  the  kth  order  time  derivatives  are 

JV+fc+l 

vL  =  C°’k+  £  (Gp cos Sf  sin. //•■)  (26) 

3  =  1 
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where 


C?’k  = 


.  Mh„o  (Pet  fM'l  ,  i  ,,WJ  „0  OM  fhh 
*  +  nk  k  \  dtk  ’  dtk  J  +  2 a[  ’  16  k  {  dtk  ’  dtk 


\ 


clk  =  i 


m.m. 


\  dtk  ’  dtk  J 


&k  =  i 


The  recursions  (27)  are  initialized  by 


£fU,U  _ 

cf  =  c? 

S i°  =  sf 


\  <9tfc  ’  <9d  / 


(27) 


(28) 


where  the  Cf,  C/,  S'/  are  the  coefficients  in  the  expansion  (25).  The  quantities  Up  Vf  in  (27) 
are  given  by  the  relations  (2.5.2-12)  with  J  =  IV  +  1  and  m  =  k ,  and  Zf  (j )  is  the  inclusion 
operator  (2.5.3-18). 


4.3  Central-Body  Gravitational  Tesserals 

For  the  central-body  gravitational  tessera!  harmonics,  the  appropriate  disturbing  function 
77  is  (2.7.1-16)  with  m  d  0.  From  the  results  in  Section  2.5.4,  we  can  construct  the  first- 
order  short-periodic  variations  r]ia.  Further  details  beyond  those  given  here  may  be  found  in 
[Proulx,  McClain,  Early,  and  Cefola,  1981]  and  [Proulx,  1981]. 

In  the  absence  of  explicit  time-dependence,  the  first-order  short-periodic  variations  are 
given  by  (2. 5. 4-4): 

OO  OO 

Via  =  E  E  lcim  cos (i A  -  mO)  +  S{m  sin(jA  -  m0)\ 

j=—oom=l  V 1  / 

The  Cf*  and  Sjm  are  given  by  (2. 5. 4-5)  in  terms  of  Fourier  coefficients  Cjm  and  S{m  of 
the  osculating  rate  functions  (2. 5. 4-1).  These  latter  coefficients  are  easily  constructed  by 
substituting  the  disturbing  function  (2.7.1-16)  with  m  / 0  into  (2.2-10).  The  needed  partial 
derivatives  of  77  have  the  same  form  as  (3.3-4).  To  illustrate,  we  give  below  the  final  formulas 
for  the  coefficients  C{m  and  S{m: 


C{m 


S 


jm 

1 


2/d 


N 


N 


A 


E  E 


s=—N  n=max(2,m,|s|) 

2«  £  £ 


(- 


A  s=— N  n=max(2,ra,|s|) 


R 


pT/mrin  j, — n—l,spvw/ptj  o  _  Zjj  f 1 

1  VnsLns1'-j  rt  \^TmsDnm  rLms'^nm 

rmirmpm  j. — n—\,spvw(p<j  pi  i  pj  c 
1  VnsLns1'-j  ri  \^rms'-ynm  ~r  rlrnsD'n 


(2) 
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(Remember  that  A  is  defined  by  (2.1.6-la).) 

This  theory  has  been  implemented  up  to  a  50x50  gravity  field  model  [Fonte,  1993]. 

4.4  Atmospheric  Drag 

For  atmospheric  drag,  the  osculating  rate  functions  Fia  are  given  by  (2.2-5)  with  perturbing 
acceleration  q  given  by  (3.4-3).  From  the  results  in  Section  2.5.1,  we  can  construct  an 
expansion  in  the  mean  longitude  A  for  the  first-order  short-periodic  variations: 

OO 

Via  =  cosiX  +  si  sinjA)  (1) 

3= 1 

where  Cj  and  Sj  are  given  by  (2.5.1-15)  in  terms  of  Fourier  coefficients  Cj  and  Si  of  the 
A-expansion  (2.5.1-11)  of  the  osculating  rate  functions: 

c‘  =  (2) 
sf  = 

Of  course,  the  limits  (Ai,  A2)  indicate  the  values  of  A  at  atmosphere  entry  and  exit,  and  we 
convert  the  integrals  over  A  into  integrals  over  F  or  L  using  (2. 5. 2-6)  or  (2.5.3-14).  Green 
[1979]  used  the  A-expansion  (1)  for  atmospheric  drag,  but  indicated  the  desirability  of  an 
expansion  in  another  geometric  variable. 

From  the  results  in  Section  2.5.2,  we  can  construct  an  alternate  expansion  in  the  eccentric 
longitude  F: 

OO 

Via  =  +  YXCi  cosJF  +  Si  sin  jF)  (3) 

3= 1 

where  Cj  and  Sf  are  given  by  (2.5.2-14)  in  terms  of  Fourier  coefficients  Cj  and  Si  of  the 
F-expansion  (2. 5. 2-1)  of  the  osculating  rate  functions: 

c{  = 
si  = 

But,  as  was  indicated  in  Section  4.1,  the  most  desirable  expansion  for  atmospheric  drag 
is  in  terms  of  the  true  longitude  L.  From  (2.5.3-23) 

K-\-  2  00 

v,a  =  c°  +  Y.  D?(L  -  \)m  +  Y.(C!  COSJL  +  S{  srnjL)  (5) 

m=  1  j= 1 


1  [F*(da.i  \ 

-Jf  I  ~qT  '  q  I  cosjF  dF 

1  rF*(doi  \  . 

-JF  I  -qT  ■  q  I  smjF  dF 
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where  Cj,  Sj,  D ™  are  given  by  (2.5.3-24)  (with  the  index  M  =  0)  in  terms  of  the  coefficients 
Cf,  V™  defined  by 

rj  1  fl2(d  at  \ 

U  =  —  ■  q  cos  iL  dL 

it  Jl  i  V  dr  ) 

5‘  =  lC{^'q)sinjLdL  <6) 

V™  =  0 

From  atmospheric  drag,  there  are  five  possible  sources  of  explicit  time-dependence: 

1.  Variations  in  the  solar  extreme  ultraviolet  (EUV)  flux. 

2.  Geomagnetic  storms. 

3.  Seasonal  latitudinal  variations  in  the  atmosphere  density. 

4.  Motion  of  the  diurnal  bulge  in  the  atmosphere. 

5.  Motion  of  atmospheric  tides  raised  by  the  Sun  and  Moon. 

Variations  in  the  solar  EUV  flux  can  cause  the  atmospheric  density  at  a  given  altitude  to  vary 
by  up  to  three  orders  of  magnitude  and  are  the  dominant  cause  of  error  in  lifetime  predictions 
for  near-Earth  satellites.  The  primary  source  of  variation  in  the  solar  EUV  flux  is  the  11-year 
sunspot  cycle,  and  it  might  be  assumed  that  variations  with  a  period  of  11  years  would  be 
too  slow  to  cause  appreciable  explicit  time-dependence  effects.  The  sunspot  cycle  is  far  from 
sinusoidal,  however.  At  the  beginning  of  each  cycle,  there  is  a  steep  rise  in  solar  EUV  flux 
which  could  conceivably  be  fast  enough  to  cause  significant  explicit  time-dependence  effects. 
In  addition,  there  is  a  secondary  variation  in  solar  EUV  flux  with  a  period  of  about  27  days 
caused  by  the  rotation  of  the  Sun,  which  brings  sunspots  alternately  into  and  out  of  view 
around  the  limb  of  the  Sun.  A  period  of  27  days  is  comparable  to  the  27.3-day  period  of  the 
Lunar  point-mass  perturbation,  which  is  known  to  have  significant  explicit  time-dependence 
effects.  Geomagnetic  storms  can  cause  large  variations  in  the  atmosphere  density  over  time 
scales  of  a  day  or  less.  The  potential  for  significant  explicit  time-dependence  effects  is  obvious. 
The  remaining  sources  of  explicit  time-dependence  effects  for  atmospheric  drag  are  thought 
to  be  too  small  or  too  slowly- varying  to  be  significant.  However,  explicit  time-dependence 
effects  from  atmospheric  drag  have  not  yet  been  studied  with  SST. 

4.5  Solar  Radiation  Pressure 

The  first-order  short-periodic  variations  r]ia  due  to  solar  radiation  pressure  are  formally  iden¬ 
tical  to  the  equations  in  Section  4.4  for  atmospheric  drag,  where  the  perturbing  acceleration 
q  is  given  by  (3.5-6).  For  solar  radiation  pressure,  the  simpler  expansions  (4.4-1)  in  the  mean 
longitude  A  or  (4.4-3)  in  the  eccentric  longitude  F  should  be  adequate.  Also,  Green  [1979] 
found  that  the  explicit  time-dependence  effects  from  solar  radiation  pressure  were  minor. 

As  we  have  seen  in  Section  3.5,  if  the  satellite  is  always  in  sunlight,  the  perturbing 
acceleration  q  can  be  derived  from  a  disturbing  function  1Z  which  is  nearly  identical  to  the 
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third-body  disturbing  function.  Hence  we  can  immediately  obtain  a  finite  Fourier  series  in 
the  eccentric  longitude  F  for  the  rjtQ  by  analogy  with  the  results  in  Section  4.2. 


5  Higher- Order  Terms 

Generally,  the  algebraic  complexity  greatly  increases  when  we  try  to  compute  higher-order 
mean  element  rates  and  short-periodic  variations.  Also,  it  is  assumed  that  higher-order 
terms  due  to  most  perturbations  are  usually  negligible.  In  this  chapter  we  report  on  those 
few  higher-order  effects  which  have  been  studied  to  date  with  SST. 

5.1  Second-Order  Aiap  and  r)iap  Due  to  Gravitational  Zonals  and 
Atmospheric  Drag 

The  second-order  mean  element  rates  Aiag  and  short-period  variations  Tjtag  clue  to  two  per¬ 
turbations  expanded  in  A  may  be  constructed  as  shown  in  Section  2.5.6.  We  can  calculate 
analytically  the  Fourier  coefficients  C{x  and  SjA  of  the  expansions  in  A  for  the  osculating 
rate  functions  Fi  i  due  to  the  central-body  gravitational  zonal  harmonics  by  substituting  the 
disturbing  function  (2.7.1-16)  with  m  —  0  into  equations  (2.2-10).  We  can  calculate  nu¬ 
merically  the  Fourier  coefficients  C{2  and  Sj2  of  the  expansions  in  A  for  the  osculating  rate 

dC^  dS^ 

functions  Fi  2  due  to  atmospheric  drag  from  (4.4-2).  The  partial  derivatives  and  -q^- 
needed  in  (2. 5. 6-5)  can  be  calculated  analytically  for  the  central-body  gravitational  zonals 
and  by  numerical  quadrature  for  atmospheric  drag. 

At  the  present  time,  the  only  terms  in  these  analytical  formulas  which  are  available  in 
the  SST  code  are: 

1.  The  J2-squared  auto-coupling  mean  element  rates  Am ,  correct  to  first  power  of  the 
eccentricity  and  with  (a, /3, 7)  replaced  by  the  explicit  formulas  (2. 1.9-3)  in  p  and  q. 
These  terms  were  constructed  with  the  MACSYMA  algebraic  utilities  of  [Zeis,  1978] 
and  [Bobick,  1981]. 

2.  The  J2-  squared  auto-coupling  short-periodic  variations  pm,  correct  to  zero  power  of 
the  eccentricity  and  with  (a,  (3, 7)  replaced  by  the  explicit  formulas  (2. 1.9-3)  in  p  and 
q.  These  terms  were  constructed  in  [Zeis,  1978]  and  corrected  in  [Green,  1979]. 

3.  The  J2/drag  cross-coupling  mean  elements  rates  Ail2,  correct  to  zero  power  of  the 
eccentricity  and  with  the  retrograde  factor  1=1  and  (ct,/3, 7)  replaced  by  the  explicit 
formulas  (2. 1.9-3)  in  p  and  q.  These  terms  were  constructed  in  [Green,  1979]. 

Green  [1979]  studied  the  second-order  effects  of  J2  and  drag  using  a  combination  of 
analytical  and  numerical  methods.  He  found  that  the  drag/  J2  cross-coupling  terms  Ai21  cause 
significant  effects  for  low  altitude  satellites.  He  found  that  Izsak’s  J2  height  correction  (1) 
applied  to  the  density  determination  in  the  formulas  (3.4-1,  3,  4)  gave  a  good  approximation 
to  the  Ai2\  terms.  The  following  expression  added  to  the  height  H  is  the  J2  short-periodic 
correction  (from  Section  4.1)  to  the  radial  distance  r: 


(1) 
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J2R 2 

4(1  —  e2)a 


sin2  i  cos  2 (/  +  a;)  +  (3  sin2  i 
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Green’s  explanation  of  why  Izsak’s  J2  height  correction  works  is  that  adding  the  J2  short 
periodics  to  the  elements  in  the  drag  osculating  rate  functions  and  then  averaging  is  equiva¬ 
lent  to  adding  the  drag/  J2  cross-coupling  terms  to  the  first-order  mean  element  rates,  if  we 
neglect  products  of  short-periodic  variations: 


<  Fi2(ai  +  riu, . .  ■ ,  «6  +  rjei)  > 


<  C2C1, . . . ,  ae)  >  +  y  < 

3 

Ai2  +  A*21 


dFi  2 
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5.2  Second-Order  r)iOLp  Cross-Coupling  Between  Secular  Gravita¬ 
tional  Zonals  and  Tesseral  Harmonics 

In  high-order  shallow  resonance  orbits,  the  tesseral  harmonics  which  contribute  the  most 
significant  short-period  motion  are  likely  to  be  those  with  degree  and  order  centered  around 
the  resonant  order.  For  such  orbits,  the  second-order  short-periodic  variations  due  to  cross¬ 
coupling  between  these  tesseral  harmonics  and  the  J2  secular  terms  may  also  be  significant. 
In  this  Section  we  outline  how  to  construct  these  critical  short  periodics.  For  further  details 
and  a  discussion  of  numerical  results  see  [Cefola,  1981]  and  [Cefola  and  Proulx,  1991]. 

For  the  present  purpose,  we  retain  in  the  expansions  (2.5.1-10)  of  the  osculating  rate 
functions  Cl  due  to  the  central-body  gravitational  zonal  harmonics  only  the  mean  element 
rates  An  given  by  (3.1-1): 

Ci  ~  An  (1) 

Furthermore,  we  completely  neglect  the  first-order  short-periodic  variations  r]n  due  to  the 
zonal  harmonics: 

Vn  ~  0  (2) 

The  osculating  rate  functions  C  2  due  to  the  tesseral  harmonics  may  be  expanded  in  the 
Fourier  series 

C2  =  y, \Cjm  cos(jA  —  md)  +  S{m  sin()A  —  m6)\  (3) 

j,m 

The  first-order  short-periodic  variations  due  to  the  tesseral  harmonics  are  then  given  by 
(2. 5. 4-4): 

r]i2  =  y  [Cf*  cos(jA  —  mO)  +  Sjm  sin(jA  —  md)]  (4) 

where,  in  the  absence  of  explicit  time-dependence,  the  C\m  and  Sjm  are  related  to  the  Fourier 
coefficients  Cjm  and  Sjm  by  (2. 5. 4-5). 


87 


To  obtain  the  second-order  cross-coupling  terms,  we  need  to  construct  the  functions 
Gn2  +  21  from  (2.3-27): 


r,  ,  ^  _  ^fdFi  i  dFi2  \  15  n 

GiY2  +  Gi2i  -  2^  (  7^hr2  +  -fa-Vrl  I  +  -^A*67ll7l2 


r=l  \  9ar 
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Substituting  (l)-(4)  into  (5)  yields 


(5) 


Gii2  +  Gi 21  «  E  \-Ci  cos(i A  -  m6)  +  s\  sin()A  - 


mt 


(B) 


where 


=  e 


5  /9A 


z")  \  9a. 


=  E 


5  '  d  A. 


V  190 

The  cross-coupling  short-periodics  are  then 
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7*12  +  '7*21  =  E  Pf  cos(jA  -  m6)  +  S'-7"  sin(jA  -  mt 


(8) 


where  the  coefficients  and  S1™  are  given  in  terms  of  C f"  and  by  the  relations 
(2. 5. 4-5).  The  partials  of  the  J2  mean  element  rates  needed  in  (7)  are  given  by  equations 

dc i™  dS 

(3.1-12).  The  partials  —~fi —  and  dla  of  the  tesseral  first-order  short-periodic  coefficients  are 

QC^m  dsjrn 

related  as  usual  to  the  partials  and  t)‘a  of  the  tesseral  osculating  rate  coefficients, 
which  may  be  obtained  by  differentiating  formulas  such  as  (4.3-2). 

Code  based  on  this  approximate  theory  has  been  developed  only  for  the  J2  secular  jm- 
daily  terms  (j  =  0  in  (4)),  with  finite  differences  used  to  obtain  the  partials  of  the  m-daily 
coefficients.  Construction  and  programming  of  a  complete  second-order  theory  for  a  double- 
averaged  perturbation  expanded  in  A,  9  has  yet  to  be  accomplished. 

6  Truncation  Algorithms 

Semianalytic  Satellite  Theory  contains  many  long  series  expansions.  Some  of  the  series 
are  infinite  and  hence  must  be  truncated,  whereas  others  are  finite  but  are  truncated  to 
reduce  the  computing  cost.  Automatic  truncation  algorithms  are  currently  used  in  the  SST 
code  for  three  of  these  series  expansions: 

(i)  The  Hansen  kernels  arc  initialized  by  the  infinite  series  representation  (2.7.3-10)  in  powers 
of  e2.  Convergence  of  the  series  has  been  investigated  by  [Proulx  and  McClain,  1988]  (see 
also  [Sabol,  1994]).  The  automatic  truncation  of  the  series  is  straightforward  and  need  not 
be  discussed  here. 


(ii)  The  expansion  of  the  mean  disturbing  function  (3.2-1)  for  a  third-body  point  mass  is 
automatically  truncated  by  a  procedure  described  by  [Long  and  Early,  1978].  In  the  Erst 
part  of  this  chapter  we  document  this  procedure. 

(iii)  The  expansion  of  the  mean  disturbing  function  (3.1-1)  for  the  central-body  zonal  har¬ 
monics  is  automatically  truncated  as  part  of  a  procedure  described  by  [Long  and  Early,  1978] 
for  the  truncation  of  the  averaged  nonresonant  central-body  disturbing  function.  However, 
nonresonant  tesserals  are  currently  excluded  from  the  averaged  equations  of  motion.  In  the 
second  part  of  this  chapter  we  propose  an  improved  automatic  truncation  algorithm  for  this 
case. 

All  other  series  expansions  in  the  SST  code  are  currently  truncated  manually  by  using 
tables  for  various  orbit  classifications  as  described  by  [Cefola,  1993] .  In  the  final  parts  of  this 
chapter  we  propose  automatic  truncation  algorithms  for  all  other  series  which  must  presently 
be  manually  truncated  by  an  experienced  user.  The  formulas  in  this  chapter  were  published 
in  [Danielson  and  Sagovac,  1995]. 

An  automatic  truncation  algorithm  removes  from  an  expansion  all  terms  whose  absolute 
values  are  less  than  a  certain  truncation  tolerance.  To  do  this,  the  algorithm  evaluates  a 
close  upper  bound  for  the  absolute  value  of  each  term  in  the  series  and  sets  the  indices  of  the 
expansion  to  include  only  those  terms  whose  upper  bounds  are  greater  than  the  truncation 
tolerance.  The  upper  bounds  are  evaluated  using  the  parameters  existing  at  the  initial  epoch 
of  the  integration  span. 

For  the  truncation  procedure  to  be  reliable,  the  upper  bounds  must  satisfy  the  following 
two  conditions  stated  in  [Long  and  Early,  1978]: 

1.  The  upper  bounds  must  be  upper  bounds  throughout  the  integration  span.  As  the  orbit 
evolves  and  the  perturbing  bodies  move,  the  absolute  value  of  each  term  in  the  series  must 
remain  less  than  or  approximately  equal  to  the  corresponding  upper  bound.  This  condition 
can  be  satisfied  by  choosing  upper  bounds  which  depend  only  on  slowly  varying  dynamic 
parameters.  For  typical  Earth  satellites  and  typical  integration  spans,  the  slowly  varying 
parameters  are  usually  the  equinoctial  orbital  elements  (a,  h,  k,p,  q),  the  distances  R3  from 
the  center  of  mass  of  the  central  body  to  the  third  bodies,  and  the  direction  cosines  ( a ,  /?,  7) 
of  the  central-body  rotation  axis.  For  eccentric  drag-perturbed  satellites  the  eccentricity  e 
may  be  a  rapidly  decreasing  parameter,  so  the  upper  bounds  should  be  increasing  functions 
of  e. 

2.  The  upper  bounds  must  eventually  monotonically  decrease  as  the  truncatable  indices  of 
the  expansion  increase.  Each  term  in  the  series  can  be  factored  into  a  product  of  constants, 
non-zero  functions  of  the  dynamic  parameters,  and  oscillating  functions  of  the  dynamic 
parameters.  The  positions  and  numbers  of  zeros  of  an  oscillating  function  vary  as  the  indices 
of  the  function  increase.  To  avoid  premature  truncation  of  the  series,  a  smooth  upper 
envelope  is  used  as  the  upper  bound  for  each  oscillating  factor. 

Our  automatic  truncation  algorithms  require  as  inputs  only  the  values  e  and  e  of  the  trun¬ 
cation  tolerances  for  the  central-body  and  third-body  gravitational  potentials,  respectively. 
In  the  current  SST  code,  it  is  suggested  to  increase  the  truncation  tolerance  e  for  distant 
eccentric  satellites  because  the  expansion  for  the  third-body  mean  gravitational  potential 
converges  slowly  for  these  satellites,  and  the  extra  accuracy  achieved  by  using  a  smaller 
tolerance  is  expensive  [Long  and  Early,  1978].  However,  implementation  of  our  proposed 
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truncation  algorithms  may  make  economical  the  use  of  a  single  relatively  small  number  for 
e  and  e  to  yield  high  accuracy. 

6.1  Third-Body  Mean  Gravitational  Potential 

The  mean  disturbing  function  due  to  the  gravitational  field  of  a  third-body  point  mass  is 
from  (3.2-1) 


S<N  N 

u  =  E  £  um  (i) 

s=0  n=max(2,s) 

where 

usn  =  §-{2  -  50s )  f^-Y  VnsK™QnsGs  (2) 

tt.3  V-K  3/ 

The  expansion  (1)  has  two  truncatable  indices:  n  giving  the  power  of  and  s  giving  the 
power  of  e.  The  purpose  of  the  truncation  algorithm  is  to  determine  the  maximum  values 
N  and  S  of  n  and  s,  respectively. 

Now  each  term  (2)  in  the  series  (1)  is  less  than  or  equal  to  the  product  of  the  absolute 
values  of  its  factors: 

usn  <  ^  \Vns\  \Kqs\  \Qns\  \Ga\ 

From  (2. 8. 1-3),  we  can  easily  obtain  the  constants  \Vns\.  From  (2.7.3-3b,5),  |iFgs(e)|  are 
clearly  positive  for  all  0  <  e  <  1.  The  functions  | (t)  I  are  bounded  for  all  —  1  <  7  <  1  by 
Qns(  1): 

( 77  -I-  S  V 
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From  (3.1-4)  and  the  relation  k2  +  h2  =  e2,  the  upper  bound  on  \GS\  is 

|GS|  <  (k2  +  h2)s/2(oi2  +  (32)s/2  =  es(a2  +  (32)s/ 2  <  es  (4) 

Multiplying  the  upper  bounds  of  all  these  factors  together,  we  finally  obtain 
Tj  <  \tj  |  _  ^3  fa)n  (n  +  s)l  fe\s 

Usn  S  \U  sn\  Bound  ~  ^  [2RJ  s,  'A0  I  [2 )  W 

The  truncation  algorithm  requires  the  calculation  of  Usn  \  Round  for  each  s  —  0, 1, . . .  and 
n  =  max(2,  s), . . .  for  n  —  s  even.  Then  N(s)  is  the  greatest  integer  for  which  Usm \ Bound  >  e, 
and  S  is  the  greatest  integer  for  which  \Usn\ Bound  >  e.  Here  e  is  the  prescribed  truncation 
tolerance  for  the  third-body  gravitational  potential.  Note  that  if  the  n  summation  in  (1) 
is  done  first  for  each  successive  s,  then  N  depends  on  s  as  indicated  above;  whereas  if  it  is 
desired  to  perform  the  s  summation  in  (1)  first  for  each  successive  n,  then  S  will  depend  on 
n. 
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This  algorithm  has  been  successful  and  is  used  in  the  current  SST  code.  However  the 
index  N  is  taken  to  be  the  maximum  N(s)  amongst  all  s,  resulting  in  the  retention  of 
negligible  terms  in  the  series  for  U .  The  upper  bound  (5)  can  be  simplified  considerably,  but 
at  a  cost  of  significantly  overestimating  the  size  of  \Usn\ Bound  dr  most  cases  [Long  and  Early, 
1978], 

6.2  Central-Body  Mean  Zonal  Harmonics 

The  mean  disturbing  function  due  to  the  gravitational  zonal  harmonics  of  the  central  body 
is  from  (3.1-3) 


S<N—2  N 
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5=0  n=s-\-  2 

where 
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The  expansion  (1)  has  two  truncatable  indices:  n  giving  the  order  of  the  geopotential 
coefficients  and  s  giving  the  power  of  e.  The  purpose  of  the  truncation  algorithm  is  to 
determine  the  maximum  values  N  and  S  of  n  and  s,  respectively. 

Now  each  term  (2)  in  the  series  (1)  is  less  than  or  equal  to  the  product  of  the  absolute 
values  of  its  factors: 
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From  (2.8. 1-3)  we  can  easily  obtain  the  constants  |V^S|.  From  (2.7.3-3a,5),  A'0"n_1’s(e) 
are  clearly  positive  for  all  0  <  e  <  1.  The  functions  |Qns(7)|  can  be  replaced  by  the  upper 
bound  \Qns(l)\ Bound,  given  in  [Danielson  and  Sagovac,  1995,  Appendix  A]: 
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From  (6.1-4)  and  the  relation  a2  +  /32  +  y2  =  1,  the  upper  bound  on  |GS|  may  now  be  taken 
to  be 


\GS\  <  es(a2  +  (32)s/ 2  =  es(l  -  72)s/2 


Note  that  since  7  is  now  a  slowly  varying  parameter,  it  need  not  be  removed  from  the 
upper  bounds  as  in  (6. 1-3,4).  Multiplying  the  upper  bounds  of  all  these  factors  together,  we 
finally  obtain 


U sn  <  |C 


2/1 


sn  Bound 


2a  J 


I  Jn 


(n 


K, 


—n—l,s 


I Q 


ns  Bound 


(1  -  72)s/V 


The  truncation  algorithm  requires  the  calculation  of  \Usn\Bound  for  each  s  —  0, 1, . . .  and 
n  =  s  +  2,  s  +  3, . . .  for  n  —  s  even.  Then  N(s )  is  the  greatest  integer  for  which  \USN\Bound  >  D 
and  S  is  the  greatest  integer  for  which  \Usn\ Bound  >  e-  Here  e  is  the  prescribed  truncation 
tolerance  for  the  central-body  gravitational  potential.  Of  course,  N  can  be  no  larger  than 
the  index  of  the  highest  available  geopotential  coefficient  J^. 
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6.3  Central-Body  Tesseral  Harmonics 

The  disturbing  function  due  to  the  gravitational  tesseral  harmonics  of  the  central  body  is 
from  (2.7.1-16) 

J2  M<N  s=S2<N  N 

R  =  E  E  E~  E  Ri™n  (1) 

j=J\  m=l  s=Si>— N  rc=max(2,m,|.s|) 

where 


Rrnm  =  Re  {£  (f  )“ ImV,T,TZK-"-'-‘Pr(C„m  -  iSnm)(Gi,  +  -  m0)}\ 

(2) 

The  expansion  (1)  has  four  truncatablc  indices:  n  giving  the  order  of  the  geopotential 
coefficients,  s,  m  giving  the  degree  of  the  geopotential  coefficients,  and  j  giving  the  frequency 
of  A.  The  purpose  of  the  truncation  algorithm  is  to  determine  the  maximum  value  N  of  n, 
the  minimum  and  maximum  values  S\  and  S 2  of  s,  the  maximum  value  M  of  m,  and  the 
minimum  and  maximum  values  J\  and  J2  of  j . 

Now  each  term  (2)  in  the  series  (1)  is  less  than  or  equal  to  the  product  of  the  absolute 
values  of  its  factors: 

\Rim,n\  <  -  (-)“ \vz\\KMr^‘Wpr\\Cnm  -  »s„m||G:„,  +  «rjj 
a  \a  )  J 

From  (2.7.1-6,13)  we  can  easily  obtain  \V™\  and  |T™  |.  The  Hansen  kernels  can  be  replaced 
by  the  upper  bound  [Danielson  and  Sagovac,  1995,  Appendix  B] 

[A'f  (e) |  <  \K™(e)\Bound  =  (1  -  e2)n+ 3/2  max  /Cf  (e) 

j  j  6—0  or  6—1  J 


where  /C”s(e)  =  |iF”s(e)|/  (1  —  e2)n+3//2.  Here  the  values  of  Kr-S  at  e  =  0  may  be  calculated 
recursively  from  (2.7.3-10,11,12,13)  or  directly  from 


0)  = 


lys  Jyl(n  -  .j  ■  /.•  -2),  ;  ,/ 


k= 0 


k\  (s  —  j  —  k)\ 


for  s  >  j 


siUn  •  k-  2),  ,  )A'  , 

'  L - 7T7Z- - —  TTi -  for  s  <  j 


k= 0 


k\  (j  —  s  —  k)\ 


where  ( a)k  are  the  Pochhammer  symbols  defined  by  (a)0  =  1  and  ( a)k  =  a(a  +  l)(a  +  2)  ■  ■  ■ 
(a  +  k  —  1) ,  and  the  values  of  /C”s  at  e  =  1  may  be  calculated  from 


— n—  2 


y(i)  =  E 

fc=0 


—n  -  2' 


0 


[l  +  (-l)fc+1  <!  2_fc_1 1 


'  k  > 
Sk-s)/ 2y 


k  <  s 
s  <  k 
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where 


'71' 


ni 


( n  —  k)\k\ 

polynomials  was  given  in  [Long  and  Early,  1978]: 


are  the  binomial  coefficients.  A  smooth  upper  bound  for  the  Jacobi 


|JT(7)I  <  I  PH  Bound  = 


\ 


[PT  (7)]2  + 


7 


d_ 

d'y 


PTh ) 


l(v  +  W  +  1  +  1) 

The  upper  bound  for  \Cnm  —  iSnm\  is  simply 

|Cnm  —  iSnm\  <  \j Cln  +  S%m 
From  (2.7.1-14)  the  upper  bound  on  | GJms  +  iH^ls |  is 

I Gjms  +  iWms\  <  ( k 2  +  h2)ls~jl/2  (a2  +  ^2)1^-7^172  =  e\s-j\^  _  ^\s-Im\/2 

Multiplying  all  of  these  upper  bounds  together,  we  finally  obtain 


(3) 


p  <£  \  p  I  _  f  |T/m| \Vm  1 1  K~n~ l,s I 

Pi'jmsn  |  J-*'jmsn  \  Bound  l  I  I  *  ns  II 1  ns  II  j  |  Bound 

CL  \  CL  / 

|fTlflo„„VQ,nH-SL(l  - 

The  truncation  algorithm  requires  the  calculation  of  | Rjmsn \ Bound  for  each  j  =  0,  ±1,  ±2, . . . 
and  m  —  1,2,...  and  s  —  j,  j  ±  1,  j  ±  2, . . .  and  n  =  max(2,  m,  |s|), ...  for  n  —  s  even.  Then 
N(j,m,s)  is  the  greatest  integer  for  which  R,jmsN \ Round  >  e,  Si(j,m)  and  ^(j, m)  are  the 
smallest  and  greatest  integers  for  which  |  Rrm,Sn  \  Bound  >  e,  M (j )  is  the  greatest  integer  for 
which  |i?jMsn|sound  >  0  and  Ji  and  J2  are  the  smallest  and  greatest  integers  for  which 
I  Rjmsn  |  Bound  >  e-  Here  again  e  is  the  truncation  tolerance  for  the  central-body  gravitational 
potential.  Of  course,  N  and  M  can  be  no  larger  than  the  indices  of  the  highest  available 
geopotential  coefficients  Cnm  and  Snm- 


6.4  Central-Body  Zonal  Harmonics  Short-Periodics 

The  short-periodic  generating  function  due  to  the  gravitational  zonal  harmonics  of  the  central 
body  is  from  (4.1-11): 


j 

5  =  U(L  -  A)  +  C°  +  J2(Cj  cos  jL  +  Sj  sin  jL)  (1) 

l=i 

Here  U  is  given  by  (6.2-1)  and 

M  2  U 

U(L  —  A)  =  ^2  — (crm  cos  mL  —  pm  sin  mL)  (2) 

m=  1  171 

Other  Fourier  coefficients  in  (1)  are  from  (4.1-12,  13,  14a) 

c°  =  -E(cVj  +  ^) 

3=1 
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Si<Ni  —  l  Ni  S2<N2~j  N2 

c‘=x E  E  Cg+Ifm-'U)  E  E  eg 

s=j  n=s+l  s=0  n=max(j+sj+l) 

Si<j-1  S5<j  N5 

+ii>iN>(j)ci3  +  ii‘iN‘U)  Y  cg+ig^-'o)  Y.  E  eg 

s=l  s=l  n=j+ 1 


^_j-Jq<2Nq  —  1 


®6<  §  —  1  min(j— 1,Nq) 

0)  Y  E  c?n 

s=j— min(j  —  1,Nq)  n=j—s 


Sr<min(j— l,Nj)—  1  min(j— l,Nj) 
.J7<2N7-)-{a\  X  '  ^  Cj7 


+ir‘"7  o)  y 


s=j/2 


n=s+l 


(3) 


where 


Cg  =  --(2  - 

aj 

eg  =  -(2  -  *,,+„)  JnH,,,+jK^-''Lg‘ 
aj 

V3  =  —.IIlujK^  ■'  u7.;;  (4) 

CL]  J 

eg  =  ij(2  - 

eg  =  Ci«  =  eg  =  -4(2  -  S„j.,)JJ,j-,K,Jn-UL i- 

aj 

For  brevity,  we  have  in  the  last  two  series  in  (3)  only  shown  the  limits  appropriate  for  j  even, 
and  we  have  not  shown  the  expansions  for  the  Fourier  coefficients  S J  in  (1).  However,  the 
maximum  values  of  the  truncatable  indices  of  the  series  we  do  not  show  are  identical  to  the 
ones  determined  by  the  procedure  outlined  here. 

The  expansion  (2)  has  the  one  truncatable  index  M  giving  the  frequency  of  L.  The 
expansions  in  (3)  have  up  to  three  truncatable  indices:  n  giving  the  order  of  the  geopotential 
coefficients,  s  giving  the  power  of  e,  and  j  giving  the  frequency  of  L.  The  purpose  of  the 
truncation  algorithm  is  to  determine  the  maximum  value  M  of  m,  the  maximum  values 
Nu...,N7  of  a,  the  maximum  values  Si, . . . ,  S7  of  s,  and  the  maximum  values  J\, . . . ,  J7  of 

j- 

Now  each  Fourier  coefficient  in  the  series  (2)  is  less  than  or  equal  to  its  absolute  value, 
which  using  (2. 5. 3-4, 5)  is 


2Uam  <  2jj7j  |  crm 

m  ~  m 


m 

m 


(1  +  mB)\b\m\Sm(k,  h)\ 


< 


m 

m 


{l+rnB)\b\m\Cm(k,h)+iSm(k,h)\ 
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(5) 


_  2|?7|  (1  +  mB )  ,  2  j^2\m/2 

m  (1  +  B)m  {  ’ 

2\U\  (1  +  myT^?)^ 
m  (1  +  a/1  -  e2)m 

The  upper  bound  of  the  Fourier  coefficient  —  2U£m  is  the  same. 

Each  of  the  Fourier  coefficients  (4)  is  less  than  or  equal  to  the  product  of  the  absolute 
values  of  its  factors,  which  using  (4.1-10a)  are 

\CH\  <  ^  | Jn\  | Hs^j\ \Vn,s-j\ Kon~hs  \Qn,s-j\ 

\c*\  <-(-)"  |/.|  |UJ+S|  K?-u  \Qn,i+,\ 

\C,3\  <  2£j  (§)' '  \Jj\  l^ojl \V,j\ K^'-°  \Qjj\ 

iq4i  <  %  (£/  u,\  ic.i-,1  im-.i  jv-1-  w«-.i 

|c£|  <  —  (-)"  |/.|  I/.J-.I  AV1'*  \QnJ~s\ 

From  (2.8. 1-3),  we  can  easily  obtain  the  constants  \Vns\.  The  functions  Kyn~l,s{e)  are 
positive,  and  the  functions  |Qn,s(7)|  may  be  replaced  by  the  upper  bound  (6.2-3).  From  the 
definitions  (2. 5. 3-5)  and  (4.1-10c),  the  upper  bound  on  HjS  is 


\Hjs\  =  |Re{ [Cj(k,  h )  +  *  Sj(k,  h)]  [S8(a,  (5)  +  *  Cs(a,  /3)]}| 

<  I  [Cj(k,  h )  +  i  Sj{k,  h)]  [S8(a,  (3)  +  i  C8(a,  (3)]  \ 

=  \Cj(k,h)  +  iSj(k,h)\  \C8(a,P)-iS8(a,P)\ 

=  (k2  +  h2yi\a 2  +  (32y ’2  =  y  (i  -  7  2y'2 

The  upper  bound  on  IjS  is  the  same.  Multiplying  the  upper  bounds  of  all  these  factors 
together,  we  finally  obtain 

Psln\  <  \Cyn\Bound  =  ^  (-)" \Jn\  Kyn-hs\Qm\Bound  (1  -  72)^  es 
\Ci2n\  <  \Cil\Bound  =  ^  (-V  I  Jn\  \vn  j+s\  Kyn~u\QnJ+s\Bound  (1  -  y2)2^  e* 

aj  \a  ) 

i cy  <  \c£\ Bmd = y  yy  w  ivy  a  -  ^ 

CLJ  \  CL  / 
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\Ci4\  <  \Ci4\Bound  =  ^  (-Y  \Jj\  \Vjj-a\  K^lQj^Bound  (1  -  72)^  es 

\Ci5n\  <  \Ci5n\Bound  =-(-)"  \Jn\  \Vn,j-s\  KoJ~1,S\QnJ-s\Bound  (1  72)^  e* 

CLJ  \  CL  / 


To  determine  the  tolerance  for  the  Fourier  coefficients  in  the  short-periodic  generator  S, 
we  note  from  (2. 5. 5-4)  that 


dS 

~d\ 


R-U 


(6) 


and  again  let  e  be  the  truncation  tolerance  for  the  central-body  gravitational  potential.  It 
follows  from  (6)  that  a  term  &  cos  jL  in  S  may  be  significant  only  if 


•  1/77  1  dL 

',|c|aX 


(7) 


Now  from  (2.5.3-14) 


d L  n - 2  (  a\2  (1  +  ecos/)2 

9A=Vl“6  \r  /  =  (1  -  e2)3/2 


From  (7)-(8) 


(1  +  e)2  _  y/l  +  e 

(1  -  e2)3/2  “  (1  -  efR 


P\  > 


(1  —  e)3/2 

jV  1  +  e 


(8) 


so  6j  is  the  appropriate  truncation  tolerance  for  the  Fourier  coefficients  in  S. 

The  truncation  algorithm  for  the  series  (2)  requires  the  calculation  of  the  upper  bound 
(5)  for  each  m  =  1,  2,  3, . . ..  Then  M  is  the  greatest  integer  for  which 


9  |77l  ^  +  e  (1  +  My/1  -  e 2) 
1  1  (1  -  e)3/2  (1  +  VT^)M 


>  e 


The  truncation  algorithm  for  the  first  series  in  (3)  requires  the  calculation  of  \C{4\Bound 
for  each  j  —  1,2,...  and  s  —  j,j  +  1, . . .  and  n  =  s  + 1,  s  +  2, . . .  for  n  —  s  even.  Then  N\  (j,  s ) 
is  the  greatest  integer  for  which  \BoUnd  >  £j,  S\  (j)  is  the  greatest  integer  for  which 
legJ  Bound  >  and  Ji  is  the  greatest  integer  for  which  \C^  \Bound  >  The  maximum 
indices  N2, ...  ,N7  and  S2, ... ,  S7  and  J2, . . . ,  J7  are  similarly  determined.  Of  course,  the 
indices  N±, . . . ,  N7  can  be  no  larger  than  the  index  of  the  highest  available  geopotential 
coefficient  Jn-  The  index  J  in  (1)  is  the  maximum  index  amongst  J±, . . . ,  J7. 

The  first-order  short-periodic  variations  are  then  given  by  equations  (4.1-18)  through 
(4.1-25)  with  the  index  N  replaced  by  -hyl. 


96 


6.5  Third-Body  Short-Periodics 


The  short-periodic  generating  function  due  to  the  gravitational  field  of  a  third-body  point 
mass  is  from  (4.2-13) 


S  —  C°  +  U(k  sinF  —  h  cos  F)  +  {C3  cos  jF  +  Sj  sin  jF) 

3= i 


(1) 


Here  U  is  given  by  (6.1-1)  and 


_n  k,  h, 
c°  =  -C1  +  -S1 
2  2 


Si<Ni  Ni  S2<N2  N2 

c‘  =  y  y  cjj+i^+'o)  y  y  eg  (2) 


s=0  n=max(2j  — l,s) 


8= 0  n=max(2j— l,s) 


where  now 


eg  =  (2  -  u  (■£■)"  u,  e-'2-> 1  TO;+1-‘ 

[sgn(j  -  5)  Ca(a,  (3)  S\j_s\(k,  h )  +  Fs(a,  /?)  C'b-_s|(/c,  /i)]  (3) 

eg  =  ^  (2  -  «  (|j)"  14,  q„»  e-0+*> 

[-Cs(a,  /?)  -%+<,(&,  /i)  +  Fs(a,  /3)  Cj+S(k,  h)] 

For  brevity,  we  have  not  shown  the  expansions  for  S J  in  (1).  However,  the  maximum  values 
of  the  truncatable  indices  of  the  series  we  do  not  show  are  identical  to  the  ones  determined 
by  the  procedure  outlined  here. 

The  expansions  in  (2)  have  three  truncatable  indices:  n  giving  the  power  of  s,  and 
j  giving  the  frequency  of  F.  The  purpose  of  the  truncation  algorithm  is  to  determine  the 
maximum  values  N,  S,  and  J  of  n,  s,  and  j,  respectively. 

Each  of  the  Fourier  coefficients  (3)  is  less  than  or  equal  to  the  product  of  the  absolute 
values  of  its  factors,  which  using  (2. 5. 3-5)  is: 

VSi  <  jg-s  (U  IK, I  K?„(X)|KMI«— 

VSi  <  jjf  (^)"  IK, I  K?„(i)|k5Ml»~- 

From  (2.8. 1-3)  and  (6.1-3)  we  can  obtain  the  constants  \Vns\  and  |(5ns(l)|.  From  (4.2-10)  we 
can  obtain  the  upper  bounds  on  |u4ls(e)|  (note  that  the  Jacobi  polynomials  P“;3(x)  >  0  for 
argument  y  >  1). 

To  determine  the  tolerance  for  the  Fourier  coefficients  in  the  short-periodic  generator  S, 
we  again  note  from  (6.4-6)  that  a  term  cos  jF  in  S  may  be  significant  only  if 

d  F 

>'c,'ax>*  <4> 


97 


where  e  is  the  tolerance  for  the  third-body  gravitational  potential.  Now  from  (4.2-18c) 


d F  a  1  +  e  cos  / 

d\  r 


< 


1  +  e 


1  —  e2 
1 


(5) 


From  (4)- (5) 


P\  > 


(1-e). 
J 


-e  =  e,- 


so  6j  is  the  appropriate  truncation  tolerance  for  the  Fourier  coefficients  in  S. 

The  truncation  algorithm  for  the  first  series  in  (2)  requires  the  calculation  of  \Ci*\ Bound  f°r 
each  j  =  1,2,...  and  s  —  0, 1, , . .  and  n  =  max(2,  j  —  1,  s), . . .  for  n  —  s  even.  Then  Ni(j,  s ) 
is  the  greatest  integer  for  which  |C^JBcmnd  >  S\  (j )  is  the  greatest  integer  for  which 
icgj  Bound  >  tj,  and  J\  is  the  greatest  integer  for  which  |C/,(J  \ Bound  >  ep  The  maximum 
indices  N2,  S2,  and  J2  are  similarly  determined.  The  index  J  is  the  maximum  index  amongst 
J\  and  J2. 

The  first-order  short-periodic  variations  are  then  obtained  from  S  by  the  procedure  out¬ 
lined  in  Section  4.2. 

6.6  Nonconservative  Short-Periodics  and  Second-Order  Expan¬ 
sions 


The  first-order  short-periodic  variations  for  a  nonconservative  perturbation  have  the  form 
(4.4-1) 


J 

Via  =  J2(Ci  COS +  Sl  Shl  J' A )  ( 1 ) 

3= 1 

The  Fourier  coefficients  in  (1)  are  determined  by  numerical  integration  of  the  osculating  rate 
functions.  The  second-order  short-periodic  variations  have  expansions  analogous  to  (1),  with 
the  Fourier  coefficients  related  to  products  of  the  osculating  rate  functions  and  first-order 
short-periodic  variations.  Although  we  have  shown  a  A-expansion  here,  it  may  be  preferable 
to  use  alternate  expansions  in  L  or  F. 

The  purpose  of  the  truncation  algorithm  is  to  determine  the  maximum  value  J  of  the 
index  j  in  (1).  We  propose  to  simply  retain  all  Fourier  coefficients  which  are  greater  (in 
absolute  value)  than  the  largest  (in  absolute  value)  Fourier  coefficient  which  has  been  dropped 
from  the  first-order  short-periodic  variations  for  the  conservative  perturbations.  This  latter 
coefficient  is  the  largest  (in  absolute  value)  of  all  the  Fourier  coefficients  in  the  short-periodic 
variations  neglected  when  applying  the  truncation  procedures  outlined  in  the  preceding  three 
sections. 

The  second-order  mean  element  rates  are  also  obtained  from  expansions  of  products  of  the 
osculating  rate  functions  and  first-order  short-periodic  variations.  We  propose  to  truncate 


these  expansions  by  retaining  all  terms  which  are  greater  (in  absolute  value)  than  the  largest 
(in  absolute  value)  term  which  has  been  dropped  from  the  first-order  mean  element  rates 
for  the  conservative  perturbations.  This  latter  coefficient  is  the  largest  (in  absolute  value)  of 
all  the  terms  in  the  mean  element  rates  neglected  when  applying  the  truncation  procedures 
outlined  in  the  first  three  sections. 


7  Numerical  Methods 

The  numerical  methods  which  are  currently  used  in  SST  are  standard.  In  this  chapter  we 
record  the  essential  mathematical  formulas.  Further  details  may  be  found  in  any  numerical 
analysis  textbook  (e.  g.,  [Ferziger,  1981]). 


7.1  Numerical  Solution  of  Kepler’s  Equation 

The  equinoctial  form  of  Kepler’s  Equation  is  (2. 1.4-2): 

A  =  F  +  h  cos  F  —  k  sin  F 


This  equation  can  be  solved  iteratively  using  Newton’s  method: 


Fo 
Fi+ 1 


/ Fi  +  h  cos  Fi  —  k  sin  Ft  —  \\ 
v  1  —  h  sin  Ft  —  k  cos  Fi  J 


for  i  —  0, 1, : 


(1) 

(2) 


7.2  Numerical  Differentiation 

We  need  to  differentiate  functions  in  order  to  obtain  the  mean  element  rates,  short-periodic 
variations,  and  partial  derivatives  for  state  estimation.  Analytical  formulas  are  preferable  if 
possible  to  obtain,  because  of  their  greater  precision.  However,  the  derivatives  of  a  function 
can  be  approximated  by  finite  difference  schemes. 

We  suppose  f(x)  is  a  smooth  function  of  x.  Then  the  central  difference  approximation 
for  the  derivative  of  f{x)  is 


dx 


(x) 


f(x  +  A)  -  f(x 
2  A 


A) 


The  error  in  this  approximation  is 


df_ 

dx 


(x) 


f(x  +  A)-f(x-A) 
2A 


A2  d3f 
6  dx 3 


(0, 


x  —  A  <  £  <  x  +  A 


(1) 

(2) 


For  example,  Green  [1979]  used  the  central  difference  approximation  (1)  to  calculate  the 
partial  derivatives  needed  for  state  estimation  (see  Section  2.6).  He  obtained  good  results 
with  a  step  size  of  A  =  10_5:r,  using  double  precision. 
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7.3  Numerical  Quadrature 

We  need  to  integrate  functions  in  order  to  obtain  the  mean  element  rates  and  short-periodic 
coefficients.  Analytical  formulas  are  preferable  if  possible  to  obtain,  because  of  their  greater 
precision.  Also,  they  are  more  computationally  efficient  since  analytical  formulas  need  to 
be  evaluated  only  once  per  mean  equations  integration  step,  whereas  numerical  integration 
requires  evaluation  at  each  abscissa  of  the  quadrature  [Long  and  McClain,  1976].  However, 
numerical  evaluation  of  integrals  of  the  type 

[  f(x)dx  (1) 

J  a 

is  mandatory  for  the  computation  of  the  mean  element  rates  and  short-periodic  coefficients 
involving  atmospheric  drag  or  solar  radiation  pressure  with  eclipsing.  Since  the  substitution 

^  =  2x  -  (a  +  b ) 
b  —  a 


transforms  the  integral  (1)  into 


(3) 


we  can  restrict  our  discussion  to  integrals  with  limits  between  —1  to  +1  without  loss  of 
generality. 

A  quadrature  formula  approximates  an  integral  by  a  weighted  sum  of  the  values  of  the 
integrand  at  points  on  the  interval  of  integration: 


ri  n 

/  f(0<%  «  $>*/(&),  -1  <  6  <  6  <  ■  •  .£n  <  1 

i=l 


(4) 


An  evaluation  of  different  quadrature  formulas  has  shown  the  Gaussian  quadrature  formulas 
to  be  generally  efficient  [Early,  1975].  The  weight  factors  wi:  for  Gaussian  quadratures  have 
been  tabulated,  and  the  abscissas  are  simply  the  zeros  of  the  Legendre  polynomial  of 
degree  n.  The  error  in  the  Gaussian  quadrature  formula  is 


/-i 


/m-  £“</(« 


i=  1 


22n+1(w!)4  rf2n/(Q 

(2n  +  l)(2n!)3  d£2n 


-1<£<1 


(5) 


A  polynomial  of  degree  2 n  —  1  is  integrated  exactly. 

The  appropriate  number  n  of  abscissas  in  the  Gaussian  quadrature  formulas  needed  for 
SST  can  vary  from  12  to  96,  depending  on  the  highest  frequency  components  contained  in 
the  function  to  be  integrated.  For  example,  Green  [1979]  found  that  if  the  first  10  pairs  of 
short-periodic  coefficients  are  to  be  retained  in  (2.5.1-13),  the  number  n  for  the  integrals 
(2.5.1-11)  must  be  at  least  48. 
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7.4  Numerical  Integration  of  Mean  Equations 


The  averaged  equations  of  motion  (1-2)  may  be  solved  with  a  Runge-Kutta  numerical  inte¬ 
gration  method.  We  consider  the  following  system  of  ordinary  differential  equations: 

dx 

S  =  f(x>‘)  (1) 

Here  x  denotes  the  column  matrix  of  mean  elements,  and  f  denotes  the  column  matrix  of 
mean  element  rates.  We  divide  the  t- axis  into  points  (f4,f2, . . .)  of  equal  width  h,  and  let 
Xj  =  x(tj).  Then  the  standard  fourth-order  Runge  Kutta  algorithm  is 


Xj+i  =  Xj  +  §(k4  +  2k2  +  2k3  +  k4)  for  i  =  1,2,3,... 
o 


where 


ki  =  f  (xi,ti) 

k2  =  f(xi  +  — k  i,ti  +  —) 

A  A 

k3  =  f(x*  +  —  k2,ti  +  —) 

k4  =  f(xj  +  Ak3,R  +  A) 
The  error  in  the  formulas  (2)  is  bounded  by 

:  fI5X, 


CA 


5  U/  -™-i 

dt 5 


(2) 


(3) 


(4) 


where  C  is  a  constant. 

Since  the  mean  element  rates  depend  only  on  slowly  varying  quantities,  step  sizes  A 
of  a  day  or  more  can  usually  be  used.  The  integrator  time  step  A  should  be  |  or  less  of 
the  minimum  period  r  of  the  oscillations  included  in  the  mean  equations  of  motion.  Some 
limitations  are  the  period  of  orbital  precession  due  to  J2  and  the  period  of  the  moon. 

Initial  values  of  the  mean  elements  a*(fi)  can  be  obtained  from  initial  values  of  the 
osculating  elements  dj(f4)  by  either  of  two  methods: 


1.  Numerically  integrate  the  VOP  equations  of  motion  over  a  time  interval  at  least  as 
long  as  the  period  of  the  largest  significant  short-periodic  effect  (usually  one  or  two 
satellite  orbits  -  see  [McClain  and  Slutsky,  1980]),  and  then  use  a  differential  correction 
procedure  to  find  the  initial  mean  elements  which  give  the  best  least-squares  fit  between 
the  SST  trajectory  and  the  Cowell  trajectory. 

2.  Use  successive  substitution  into  the  near-identity  transformation  (1-1)  until  a  specified 
agreement  is  reached: 

a°(U)  =  diiti) 

i)  =  ai(ti)  -  77i[aJ(ti),  •  •  ■  ,a£(ti),ti]  for  k  =  0, 1,  2, . . . 

This  method  is  faster  than  method  1,  but  may  require  the  inclusion  of  a  comprehensive 
set  of  short-periodic  variations  to  avoid  a  large  bias  in  the  initial  mean  elements. 
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It  should  be  pointed  out  that  the  time  averages  of  the  osculating  elements  over  some  time 
interval  are  generally  not  a  good  approximation  to  the  mean  elements  [Early,  1986]. 

7.5  Interpolation 

Since  the  mean  elements  and  short-periodic  coefficients  are  slowly  varying,  their  values  at 
desired  times  not  coinciding  with  the  mean  equation  step  times  can  be  computed  by  relatively 
low  order  interpolation  formulas. 

First,  suppose  that  at  distinct  times  (fi, . . .  ,tn)  we  know  the  values  [f(t i), . . . ,  f(tn )]  of 
a  smooth  function  /(f).  In  Lagrange  interpolation  we  approximate  /(f)  by  a  polynomial  of 
degree  n  —  1  passing  through  the  known  values: 

n 

f(t)  ~  /&)£*(*) 

i=  1 

Here 

L  (f)  =  (f  -  fi)  •  •  •  (f  -  -  ti+i)  •  •  •  (f  -  tn) 

( U  -  ti)  ■  ■  ■  (ti  -  -  ti+ 1)  •  •  •  (h  -  tn) 

Note  that 

Li(tj)  =  Sij  (3) 

The  error  in  the  Lagrange  interpolation  formula  is 

f  1 

Lagrange  interpolation  is  currently  used  to  interpolate  the  short-periodic  coefficients,  the 
velocity  vector,  and  the  partial  derivatives  needed  for  differential  correction.  An  adequate 
order  n  —  1  has  been  found  to  be  3  (4  interpolator  points)  [Taylor,  1978]. 

Next,  suppose  that  at  distinct  times  (ti, . . , ,  tn)  we  know  both  the  values  [/(ti), . . . ,  f(tn)] 
and  the  derivatives  [/(ti), . . . ,  f(tn)]  of  a  smooth  function  /(f).  In  Hermite  interpolation  we 
approximate  /(f)  by  a  polynomial  of  degree  2 n  —  1  passing  through  the  known  values  and 
derivatives: 


n 

fit )*  £{[i  -  2 (f  -  ti)Li(ti)][Li(t)]2f(ti)  +  (f  -  ti)[Li(t)}2f(ti)}  (5) 

1=1 


Here  again  Lj(f)  are  the  Lagrange  basis  functions  (2).  The  error  in  the  Hermite  interpolation 
formula  is 


[(t  —  ti)  •  •  •  (t  —  f«)]2  d2nf 

(2n)!  r//-’"  U  j 


fi  <  /  <  tn 


(6) 


Hermite  interpolation  is  currently  used  to  interpolate  the  mean  elements  and  the  position 
vector.  An  adequate  order  2 n  —  1  has  been  found  to  be  5  (3  interpolator  points). 
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